applyrule([

cos(u::anything)*cos(v::anything)-sin(u::anything)*sin(v::anything)=cos(u+v),

cos(u::anything)*sin (v::anything)+sin(u::anything)*cos(v::anything)=sin(u+v),

sin(u::anything)*sin(v::anything)-cos(u::anything)*cos(v::anything)=-cos(u+v),

-sin(v::anything)*cos(u::anything)-sin(u::anything)*cos(v::anything)=-sin(u+v)], simplify(-(sin(psi[1](t))*cos(theta[1](t))*cos(gamma[1](t))+sin(psi[1](t))*sin(theta[1](t))*sin(gamma[1](t))-cos(theta[1](t))*cos(psi[1](t))*sin(gamma[1](t))+cos(psi[1](t))*sin(theta[1](t))*cos(gamma[1](t)))*cos(beta[1](t)), size)):

(trigsubs(2*combine(%))[])/2;

                         cos(beta[1](t))*sin(-theta[1](t)-psi[1](t)+gamma[1](t))

There are infinitely many solutions. They all lie on the circle with center  -I  and radius 1. The equation of the circle in the complex plane will:

x:=a+I*b:

evalc(abs(I+x) = 1);

                          (a^2+(b+1)^2)^(1/2) = 1

The real solution is only one  x=0

Addition: This result is clear without any calculations because the equation can be rewritten as  abs(x-(-I))=1  and we know that  abs(z1-z2)  is the distance  between  z1  and  z2  in the complex plane.

Perhaps you work in a classic interface, in which the filling for non-convex polygons is not working properly. In the standard worksheet everything correctly:

L := [[[0, 0, 0], [0, 1, 0]], [[x, x^2+1, 0], x = 0 .. 2], [[2, 5, 0], [2, 2, 0]], [[x, x, 0], x = 2 .. 0]]:

plots[display](Picture(L, color = yellow), axes = normal, scaling = constrained);

Use  solve  instead of  fsolve:

restart;

solve({6.268 = a*c+0.710, 10.049 = a*(c+2.5*10^(-5)*(0.004/(0.025)))+0.710}, {a, c});

                                       {a = 945250., c =0 .5879925946e-5}

We can easily find maximum symbolically, if to convert the sum into the product and use that the maximum and minimum of  cos(a*t+b) , a<>0,  are the same as  cos(t) :

restart;

f:= cos(2*t/m) + cos(2*(t+5)/m):

simplify(applyrule(cos(alpha::anything)+cos(beta::anything)=2*cos((alpha+beta)/2)*cos((alpha-beta)/2), f));

max(eval(%, op(2,%)=-1), eval(%, op(2,%)=1));

restart;

eq1:=-sin(alpha0(t))*cos(beta0(t)) = -sin(alpha[1](t))*cos(beta[1](t));

eq2:=cos(alpha0(t))*cos(beta0(t)) = cos(alpha[1](t))*cos(beta[1](t));

simplify(eq1/eq2, trig);

convert(%, tan)*(-1);

This variant works immediately:

restart;

question1 := b(t)*(diff(c(t), t))*(diff(a(t), t))+a(t)*(diff(b(t), t)):

eval(question1, b(t)=0);

                                                0

See help on  Student[MultivariateCalculus][Jacobian]  command.

In your example, I have changed the variables so that they to be the names because

type(beta(t), name);

whattype(beta(t));

                                                 false

                                               function

Solution: 

eq1:=l1*cos(theta)+l2*sin(beta)-x=0:

eq2:=l1*sin(theta)-l2*cos(beta)=0:

Student[MultivariateCalculus][Jacobian](lhs~([eq1,eq2]), [beta,x,theta]);

T:=[alpha0(t)=alpha, beta0(t)=beta,psi[1](t)=psi, theta[1](t)=theta, gamma0(t)=gamma]:

A:=subs(T, eq_liaison):

add(applyrule([cos(u::anything)*cos(v::anything)-sin(u::anything)*sin(v::anything)=cos(u+v), cos(u::anything)*sin(v::anything)+sin(u::anything)*cos(v::anything)=sin(u+v), sin(u::anything)*sin(v::anything)-cos(u::anything)*cos(v::anything)=-cos(u+v), -sin(v::anything)*cos(u::anything)-sin(u::anything)*cos(v::anything)=-sin(u+v)], simplify(coeff(lhs(A),i), size))*i, i={xb[1],yb[1],zb[1],l2[1],l3[1]}):

B:=subs(map(t->rhs(t)=lhs(t),T), %);  # Final result

 Addition:

length(B);

                                          753

length(eq_liaison);

                                         1686

@Kimkimkim  Here is Konigsberg graph with 2 fictitious vertices  a  and  b , created as Carl has offered:

with(GraphTheory):

G:=Graph(Trail(1,2,3,4,1,3),Trail(1,a,2),Trail(1,b,4));

DrawGraph(G, style=planar);

@Markiyan Hirnyk   This example can be solved much shorter and with greater accuracy:

evalf(Int(exp(-a^3-b), [b=1-a..1+a, a=0..10]));

                     0.3707213163

You can use a caption instead of legend.

Example:

plot3d([(1/2)*x+(1/2)*y, 1+(1/2)*x+(1/2)*y], x = 0 .. 3, y = 0 .. 3, color = [yellow, blue], style = surface, axes = normal, orientation = [-30, 70], caption = "First plane is yellow, second plane is blue", captionfont = [Times, Roman, 14]);

Another variant - you can write the equations of the planes directly on the plot:

restart;

A := plot3d([(1/2)*x+(1/2)*y, 1+(1/2)*x+(1/2)*y], x = 0 .. 3, y = 0 .. 3, color = [green, blue], style = surface, axes = normal, orientation = [-30, 70]):

B := plots[textplot3d]([[3.3, 3.3, 2.5, z = (1/2)*x+(1/2)*y, color = green], [3.3, 3.3, 3.5, z = 1+(1/2)*x+(1/2)*y, color = blue]], font = [Times, Roman, 12]):

plots[display](A, B);

Edited.

Maple computes the integral right. But if you do not quite trust to Maple, you can do it step by step with a check these steps.

I dropped the constant factors as they do not affect the convergence of the integral. Also, we assume that  t> 0. The integrand in the interval  s=0 .. t  has only one singularity if  s=t. Therefore, we can consider the original integral as the limit of ordinary definite integrals on the range  s=0 .. p  when the variable  p  tends to  t .

A:=int(s^(5/2)/(t^(1/2)-s^(1/2)), s=0..p)  assuming t>0, p>0, p<t;

limit(A, p=t) assuming t>0;

Validation of calculation of the integral:

int(s^(5/2)/(t^(1/2)-s^(1/2)), s)  assuming t>0, s>0, s<t;

simplify(diff(%, s))  assuming t>0, s>0, s<t;

Without  solve - a line as an intersection of 2 planes:

L := [(x-2)*(1/3), (y-1)*(1/4), (z-3)*(1/3)]:

plots[intersectplot](L[1]-L[2], L[2]-L[3], x = -4 .. 4, y = -4 .. 4, z = -4 .. 4, linestyle = 1, thickness = 3, axes = normal, orientation = [60, 75], view = [-1 .. 3, -4 .. 2, -1 .. 4]);

Addition:  These planes  L[1]-L[2]=0  and  L[2]-L[3]=0  are the projecting planes on the coordinate planes x0y  and  y0z, respectively (the first one is perpendicular to x0y  and the second one is perpendicular to y0z). The equation  L[1]-L[2]=0  in the plane  x0y only, is the projection the original line into the plane  x0y  and  L[2]-L[3]=0  in the plane  y0z only, is the projection the original line into the plane y0z .

For  numeric  option, the parameters must be specified.

This version works for a small range:

restart;

dsys:={diff(u(t),t)=-(N1+M*cos(2*I*omega*t))*u(t)-1+v(t)*exp(-2*I*omega*t)+w(t)*exp(2*I*omega*t), diff(v(t),t)=-(N1+I*Delta-2*M*exp(2*I*omega*t))*v(t)-(N1+u(t))*exp(2*I*omega*t)-2*M, diff(w(t),t)=-(N1-I*Delta-2*M*exp(-2*I*omega*t))*w(t)-(N1+u(t))*exp(-2*I*omega*t)-2*M}:

N1,M,Delta,omega:=1,2,3,4:

res:=dsolve(dsys union {u(0)=-1,v(0)=0,w(0)=0},numeric,output=listprocedure, maxfun=500000);

plots[odeplot](res,[t,Re(w(t))],t=0..5,axes=boxed,color=black,linestyle=1,tickmarks=[3, 4],thickness=2);

 Addition:  Setting  maxfun=0  is dangerous if a singularity occurs, as in this case, the program can run forever. I think in your example of just such a situation.