A more detailed analysis shows that the system has 2 real solutions. The first solution, if  T__1s(0) = 62.68925412  (which Joe found) and the second one if  T__1s(0) = -62.68410750 :

restart; 
Sys := {Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (0.1375e-1*(T__1(t)-T__1s(t)))*((T__1(t)+T__1s(t))*(1/2)), diff(Q(t), t) = (0.1375e-1*(T__2s(t)-T__2(t)))*((T__2s(t)+T__2(t))*(1/2)), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2)}: 
Sys1 := subs({Q(t) = 0, T__1(t) = T1, T__1s(t) = T1s, T__2(t) = T2, T__2s(t) = T2s, diff(Q(t), t) = DQ}, Sys): 
R := eliminate(Sys1, {DQ, Q, T1, T2}); 
solve(R[2]);  # We see that specification T__1s(0)  or  T__2s(0)  is mandatory

E_T := (2/mu-2/r)*exp(-r/mu)*Pi^2;
Q1:=content(numer(select(t->type(t,`+`), E_T)));
Q2:=select(t->type(t,realcons), E_T);
Q1*Q2*combine(E_T/Q1/Q2);
                               

A simple procedure Triangle for the given side lengths of the triangle finds all its angles and draws the triangle together with the labels of the vertices and side lengths:

restart;
Triangle:=proc(a, b, c)
local AngleA, AngleB, AngleC, A, B, C, T, Tr; 
uses geometry, plots;
AngleA:=arccos((b^2+c^2-a^2)/2/b/c);
AngleB:=arccos((a^2+c^2-b^2)/2/a/c);
AngleC:=arccos((a^2+b^2-c^2)/2/a/b);
point(A,0,0);
point(B,c,0);
point(C,b*cos(AngleA),b*sin(AngleA));
print(`The angles A, B, C are:`);
print(AngleA, AngleB, AngleC);
T:=textplot([[c/2,0,c, align=below], [(c+b*cos(AngleA))/2+a/30,b*sin(AngleA)/2,a, align=right], [b*cos(AngleA)/2-a/30,b*sin(AngleA)/2,b, align=left]]);
display(draw([triangle(Tr,[A,B,C]), A, B, C], font=[times,roman,18],labels=[x,y], axes=none, printtext = true), T);
end proc:
 

Example of use:

Triangle(5, 7, 3);

Addition:  A similar procedure can easy be written, if  two sides of a triangle and an angle between them are known (as Carl did), or a side and  two angles adjacent to it, etc.

Two ways for solving the the problem.

First way (symbolic solution):

restart;
ode:=diff(y(t), t$3)+3*(diff(y(t), t$2))+4*(diff(y(t), t))+12*y(t) = 0: 
dsolve({ode, y(0)=3, D(y)(0)=0, (D@@2)(y)(0)=0});
assign(%);
plot(y(t), t=-1..5, -5..5);

Second way:

restart;
ode:=diff(y(t), t$3)+3*(diff(y(t), t$2))+4*(diff(y(t), t))+12*y(t) = 0: 
ivp := [(D@@2)(y)(0) = 0, D(y)(0) = 0, y(0) = 3]:
DEtools[DEplot](ode, y(t), t = -1 .. 5, y(t) = -5 .. 5, [ivp]);

The second method is used for the numerical solution with additional capabilities.

Maple does not have built-in commands for working with Fourier series. But you can download the app for this purpose here

Carl, the question was about the rotation of a loop of the curve . For example, a circle,  wholly located above the axis  Ox , rotates around the axis Ox. As a result, we get a torus. To get the volume of the torus, we must from body volume by rotation only the upper semicircle take away the body volume  by rotation only the lower semicircle ( r is the radius of the circle centered at the point  (0, R) ,  r<R ):

Sol:=solve(x^2+(y-R)^2=r^2, y);
f1:=unapply(Sol[1], x); f2:=unapply(Sol[2], x);
Int(Pi*f1(x)^2, x= -r..r) - Int(Pi*f2(x)^2, x= -r..r);
value(%) assuming r>0, R>0;

In  Student package:

Student[VectorCalculus][SurfaceInt](1, [x, y, z] = Surface(<r*cos(t), r*sin(t), r^2>, t = 0 .. 2*Pi, r = 0 .. 3));

By  int  command in polar coordinates:

F:=x^2+y^2:
int(r*eval(sqrt(1+diff(F, x)^2+diff(F, y)^2), {x=r*cos(t), y=r*sin(t)}), [r=0..3, t=0..2*Pi]);
expand(%);

These variants will suit you?

diff(f(x,y), x),  diff(f(x,y), y),  diff(f(x,y), x,y),  diff(f(x,y), x$2,y$3);

As you want:

`&PartialD;`*f(x,y)/cat(`&PartialD;`,x);
`&PartialD;`^2*f(x,y)/(cat(`&PartialD;`,y)*cat(`&PartialD;`,x));

Addition.

The simple procedure called  PartialD  solves your problem for an arbitrary number of variables:

restart;
PartialD:=proc(f, L)
local n, L1, L2;
uses ListTools;
L1:=Reverse([args[2..-1]]);
n:=nops(L1);
L2:=map(t->cat(`&PartialD;`,t),L1);
`&PartialD;`^n*f/convert(L2,`*`);
end proc:

Examples of use:

PartialD(f(x,y,z), x$2,y$3,z);

The arguments for the function, you don't need to write:

PartialD(H, x$2,y$3,z);

Edited.

error cannot be a name. A workaround:

restart;
`error`:=14/45*(D@@5)(y)(0)*h^5+7/10*(D@@6)(y)(0)*h^6;
remove(has, `error`, h^6);
                                      

For the plotting you must specify all the parameters. 

Example of the plotting:

restart;
schro := -(diff(psi(x), x, x))+(2*a*b*x^4+a^2*x^6+(b^2-a*(2*p+3))*x^2-(2*p+1)*b)*psi(x) = 0:
dsolve(schro);
assign(%):
plot(eval(psi(x), {a=1,b=1,p=1,_C1=2,_C2=1}), x=0..2);

multiply the latter matrix by 1/a^(3/4) :

A:=<a,b,c,d; e,f,g,h; i,j,k,l; m,n,o,p>;
B:=<a,b,c,d; 0,-e*b+a*f,-e*c+a*g,-e*d+a*h; 0,a*j-i*b,a*k-i*c,a*l-i*d; 0,-m*b+a*n,a*o-m*c,a*p-m*d>; 
D1:=LinearAlgebra[Determinant](A);
D2:=LinearAlgebra[Determinant](1/a^(3/4)*B);
D1-D2;

Download M.mw

plots[inequal]  command works for inequalities with 2 variables only. In your example 4 equalities with 3 variables are given. You can display them by plots[implicitplot3d]  command:

plots[implicitplot3d]([x-y-z=3, 2*x-y+2*z=5, x-y+z=1, 3*x-y-5*z=5], x=-7..7,y=-7..7,z=-7..7, color=[red,blue,green,yellow], style=surface, axes=normal,view=[-8.4..8.4,-8.4..8.4,-8.4..8.4]);

plots[inequal]  command works for inequalities with 2 variables. In your example 4 equalities with 3 variables are given. You can display them by   plots[implicitplot3d]  command:

plots[implicitplot3d]([x-y-z=3, 2*x-y+2*z=5, x-y+z=1, 3*x-y-5*z=5], x=-7..7,y=-7..7,z=-7..7, color=[red,blue,green,yellow], style=surface, axes=normal,view=[-8.4..8.4,-8.4..8.4,-8.4..8.4]);

Inequality  A>B  is equivalent to the inequality  A-B>0  and it will be true if

Optimization[Minimize](A-B, {constraints})[1] > 0

Examples:

A := 2*x^2+(y-2)^2-1: 
B := x^2:
Optimization[Minimize](A-B);  # On the whole plane (x,y) the inequality  A>B  does not hold 
Optimization[Minimize](A-B, {2*y <= x});  # At a given halfplane   2*y <= x  the inequality  A>B holds
is(%[1]>0);

Another way:

type(390/7, integer);

                                   false