m:=1: w[d]:=2: w[n]:=3: zeta:=4:  # I took some values of parameters
g:=t->1/m/w[d]*exp(-zeta*w[n]*t)*sin(w[d]*t);
plot(g(t), t=0..8);
f:=t->piecewise(t<0,0,t<1,100,t<2,200,t<3,5);
plot(f(t), t=0..8);
x:=t->Int(f(tau)*g(t-tau), tau=0..t):
plot(x, 0..8);

b:=sort(x^2+1, order=plex(x), ascending):
a:=sqrt(b);

or after the calculation

restart;
a:=sqrt(1+x^2);
a:=applyop(sort, 1, a, order=plex(x), ascending);
 

In  plots[fieldplot]  command the first argument should be vector field (usually the list of 2 expressions in x and y variables). Here is the example of plotting the gradient of your function:

restart;
with(plots):
A:=1.0:   z := A*x*y: 
fieldplot([diff(z,x), diff(z,y)], x=0..1, y=0..1);

Because your points are uniformly in the first coordinate, then just increase step in your data. In the example in the second plot the points along the sinusoid are 5 times less frequently:

restart;
r:=rand(-0.1..0.1):
A:=plot(sin(x), x=0..2*Pi, color=blue):
X:=[seq(2*Pi/100*k, k=0..100)]:
Y:=[seq(sin(2*Pi/100*k)+r(), k=0..100)]:
B:=plot(X, Y, style=point, symbol=solidcircle, color=red, symbolsize=7):
C:=plot([seq(X[i], i=1..101,5)], [seq(Y[i], i=1..101,5)], style=point, symbol=solidcircle, color=red, symbolsize=7):
plots[display](A, B, scaling=constrained);
plots[display](A, C, scaling=constrained);

Your calculation is erroneous, because at a certain interval the given function takes  negative values (should be  r(theta)>=0 ).

Here is a correct solution:

r:=theta->3*cos(theta)-2*sin(theta):
solve(r(theta)>=0);  # The domain of the function r
R:=op(1,%)..op(2,%);
plot(r(theta), theta=R, coords=polar);
Area:=int(1/2*r(theta)^2, theta=R);
evalf(Area);
Length:=int(sqrt(r(theta)^2+diff(r(theta),theta)^2), theta=R);
evalf(Length);

Maple is right!  

arctan(y, x)  command returns the angle which forms the radius vector of a point  with coordinates (x, y) with the positive direction of x-axis.

Always    -Pi<arctan(y, x)<=Pi

If  x>0  then  arctan(y, x)=arctan(y/x)
If  x=0  and  y>0  then  arctan(y, x)=Pi/2

If  x=0  and  y<0  then  arctan(y, x)=-Pi/2

If  x<0  and  y>=0 then  arctan(y, x)=Pi+arctan(y/x)

If  x<0  and  y<0 then  arctan(y, x)=-Pi+arctan(y/x)

Examples of use:

arctan(1, 1),  arctan(1, -1),  arctan(-1, 1),  arctan(-1, -1),  arctan(0, -1);

                            Pi/4,  3*Pi/4,  -Pi/4,  -3*Pi/4,  Pi

Edited.

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);