u1 := vector([-1,2,1,1]); u2 := vector([2,-1,-1,1]); u3 := vector([2,1,-1,3); v1 := evalm(u1/sqrt(7)); v2 := evalm( u2 - innerprod(u2,v1)*v1 ); v2 := evalm(v2/sqrt(sum(v2[i]^2,i=1..4))); sol := innerprod(u3,v1)*v1 + innerprod(u3,v2)*v2; sol := evalm(sol); A := augment(u1,u2); SOL := evalm( A &* leastsqrs(A,u3) ); B := matrix(3,3,[1,0,0,0,1,2,0,2,1]); eigensys := eigenvects(B); G := diag(-1,3,1); P := augment([0,-1,1],[0,1,1],[1,0,0]); evalm(P &* G &* inverse(P)); ExpB := evalm(P &* diag(exp(-1),exp(3),exp(1)) &* inverse(P)); evalf(");