See  Statistics[Regression]  help page. 

I could not find the options for the exact solution. But the exact solution can be easily obtained by the simple formulas:

with(Student[LinearAlgebra]):

a:=<-2|3|2>;  b:=<7|-3|-4>;

Pr:=(a.b)/(Norm(b)^2)*b;   # exact projection

evalf(Pr);   # approximate projection

c:=a-Pr;  # exact orthogonal complement

evalf(c);   # approximate orthogonal complement

Norm(c);  # exact norm of orthogonal complement

evalf(%);   # approximate norm of orthogonal complement

 Addition:  See http://en.wikipedia.org/wiki/Vector_projection

Write as follows:

g := (1/8)*x*sin(13/x)+arcsin(5*x^2/(8*x^2+1));

G := unapply(g, x);

H := unapply(G(G(x)), x);

D(H)(Pi/2);

In calculating the derivative of a function at a point convenient to use the differential operator  D  command instead of  diff  command. Here are 2 equivalent variants with  D  command:

f:=x->(x^12-x*sin(x^11))/x^34+exp(sqrt(x+4))*ln(abs(cos(x)^5-6));

D[1,1](f)(5);

(D@@2)(f)(5);

Compare with the result of work  diff  and  eval  command:

eval(diff(f(x),x$2), x=5);

Example:

Ugen := .1049253920*phi(x)^2+.2490325160*psi(x)^2+0.7218836157e-1*eta(x)^2+0.4163942054e-1*(diff(psi(x), x))^2+0.3590610475e-1*(diff(phi(x), x))^2+0.1916547983e-2*psi(x)*(diff(phi(x), x, x))+0.5733315777e-2*(diff(psi(x), x))*(diff(phi(x), x))-0.3273703041e-1*phi(x)*psi(x)-0.3273703035e-1*psi(x)*eta(x)+0.4980952543e-2*(diff(phi(x), x, x))^2-0.7191876251e-1*phi(x)*eta(x)+0.1153177175e-2*eta(x)*(diff(phi(x), x, x))+0.1073591707e-1*phi(x)*(diff(phi(x), x, x)):

coeff(subs(diff(phi(x),x)*diff(psi(x),x)=t, Ugen),  t);

                                 .5733315777e-2

Maple can not calculate your integral symbolically, but it can do this  numerically for any  n :

f := 2*sin(Pi-Pi*exp(-t)):

A[0] := (1/3)*(int(f, t = 0..3));

A := n->(2/3)*(int(f*cos((2/3)*Pi*n*t), t = 0..3));

seq(evalf(A(n)), n = 1..10);

Addition:  If you want to expand your function into Fourier series  only with cosines, it means that you continue this function from segment [0,3] as even function. I corrected some errors in the formulas. Plotted the original function and series expansion with the first 8 terms:

restart;

f := 2*sin(Pi-Pi*exp(-abs(t))): # by abs the function expanded as even function

A := n->(2/3)*evalf(Int(f*cos((1/3)*Pi*n*t), t = 0..3));

A(0)/2+add(A(n)*cos((1/3)*Pi*n*t), n=1..7);

plot([f, %], t=-3..3, color=[red,blue], thickness=2);

Probably the questioner meant getting real solutions depending on a real parameter m. But Maple does not solve the problem:

restart;

assume(m::realcons);

RealDomain[solve](x^2-(2*(m+1))*x+m^2-2*m+m^2 = 0, x, parametric = full);

                     [{x = m+1+sqrt(-m^2+4*m+1)}, {x = m+1-sqrt(-m^2+4*m+1)}]

Mathematica solves the problem correctly:

Reduce[x^2 - 2*(m + 1)*x + m^2 - 2*m + m^2 == 0, x,
Reals] 

(m == 2 - Sqrt[5] &&
x == 3 - Sqrt[5] - Sqrt[
1 + 4 (2 - Sqrt[5]) - (2 - Sqrt[5])^2]) || (2 - Sqrt[5] < m <
2 + Sqrt[5] && (x == 1 + m - Sqrt[1 + 4 m - m^2] ||
x == 1 + m + Sqrt[1 + 4 m - m^2])) || (m == 2 + Sqrt[5] &&
x == 3 + Sqrt[5] - Sqrt[1 + 4 (2 + Sqrt[5]) - (2 + Sqrt[5])^2])

Explanation:  &&  is logical  and,  ||  is logical  or

S:="110100100011001001010110001000":

L:=[1]:

for i from 2 to length(S) do

if S[i]=S[i-1] then L:=subsop(nops(L)=L[nops(L)]+1, L) else L:=[op(L), 1] fi:

od:

L; 

                               [2, 1, 1, 2, 1, 3, 2, 2, 1, 2, 1, 1, 1, 1, 2, 3, 1, 3]

Addition  -  another variant for long strings:

restart;

S:=StringTools[Random](10^5, 'binary');

L[1]:=1: k:=1:

for i from 2 to length(S) do

if S[i]=S[i-1] then L[k]:=L[k]+1 else k:=k+1: L[k]:=1 fi:

od:

convert(L,list);

fsolve command does not solve the equation with unknown parameters. Use solve command:

eq3 := k*y^2+x^2 = 4:

eq4 := -h*x^2+a*y = 0:

solve({eq3, eq4}, {x, y});

allvalues(%);

StringTools[Random](30, 'binary');

StringTools[CharacterFrequencies](%);

                     "110100100011001001010110001000"

                                  "0" = 18, "1" = 12

for n while ithprime(n)<29 do

end do:

n;

                       10

You can immediately get a decimal answer, if you use the  fsolve command, or solve command but at least one of the coefficients of your equation has float type:

A:=plot(sqrt(x), x=-1..5, thickness=3, color=blue):

B:=plot(sqrt(x), x=1..4, filled=true, color=yellow, view=[-1..5, -1..3]):

plots[display](A,B, scaling=constrained);

solve(z^3=lambda^3, z);

is(-1/2+1/2*I*3^(1/2)=exp(2*I*Pi/3));

 Addition: You can also write it in exponential form or by  j :

Sol:=solve(z^3=lambda^3, z);

z1, z2:=op(1,Sol[2]), op(1,Sol[3]);

op(subs({z1=abs(z1)*exp(``(Pi*argument(z1))),z2=abs(z2)*exp(``(Pi*argument(z2)))}, [Sol]));

op(subs({z1=j,z2=j^2}, [Sol]));

Enough to check that any vector of the subspace  V  can be uniquely decomposed into the specified basis:

V, e1, e2, e3:=<2*a+b+c, 3*a+b+2*c, 2*a+c, b+7*c>, <2,3,2,0>, <1,1,0,1>, <0,1,1,6>;

Equate(V, x*e1+y*e2+z*e3);  #  The system of linear equation

solve(%, {x,y,z});