Error reason is that the syntax of  VectorCalculus:-int  is different from the  int  command  from  Maple kernel, which you have used. Therefore, do not download the whole package  VectorCalculus, and use the long form of command call as  VectorCalculus:-command

I corrected only the first error in the calculation of the volume. The remaining errors you can easily fix yourself.

See the corrected file.

error_integral1.mw

Slightly reduce the right end of  x - range:

f:=x->piecewise(x>-Pi and x<=Pi, 3*x):
plot(f(x), x=-Pi..Pi, discont=true):
a:=-Pi:  b:=Pi:
p:=b-a:
fp:=x->f(x-floor((x-a)/p)*p):
plot(fp(x), x=-6*Pi..7*Pi-0.001, discont=true, tickmarks=[piticks,piticks]);

Example:

A:=plot3d([cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)], phi=-arccos(0.9)-0.01 .. Pi/2-arccos(0.9), theta=arccos(0.9) .. `if`(phi>=-arccos(0.9) and phi<=arccos(0.9), arcsin(0.9/cos(phi)), Pi/2), style=surface, color=khaki, scaling=constrained, numpoints=25000):
B:=plots[display](seq(plottools[rotate](plots[display](A), Pi/2*k, [[0,0,0], [0,0,1]]), k=0..3)):
C:=plottools[reflect](B, [[0,0,0], [1,0,0], [0,1,0]]):
Cube:=plottools[cuboid]([-0.9,-0.9,-0.9], [0.9,0.9,0.9], color=pink, transparency=0.8):
plots[display](B, C, Cube, axes=none);

Edit.
 

For example,  if  i=4  then

`or`(seq(R[i]=R[k], k=1..i-1))   means

R[4]=R[1]  or  R[4]=R[2]  or  R[4]=R[3]

In order to more clearly show that your function is unbounded in the neighborhood of the origin (and the top and bottom), it is helpful to reduce the ranges for x and y axes, and to increase grid and to equate the scales along the axes.

plot3d(2*x/(x^2+y^2), x = -5 .. 5, y = -5 .. 5, style=surface, color=khaki, view=-10..10, grid=[500,500], scaling=constrained, axes=normal);

If I understand the problem, the walk is over, if we find ourselves at a point where we already were in this walk. I completely rewrote your procedure  Randwalk3 . Procedure  RandWalk  describes only one random walk. To specify several walks, we simply repeat this procedure.

RandWalk:=proc(N)
local P, i, R3, r;
P[1]:=[0,0]; P[2]:=[1,0];
R3:=rand(1..3);    # 1: go straight on; 2: turn right; 3: turn left
for i from 3 to N do
r:=R3();
      if r=1 then P[i]:=[P[i-1][1], P[i-1][2]+1]      # go straight on
      elif r=2 then P[i]:=[P[i-1][1]+1, P[i-1][2]]   # turn right
      else P[i]:=[P[i-1][1]-1, P[i-1][2]]                 # turn left
      end if;
      if `or`(seq(P[i]=P[k], k=1..i-1))  then  break end if; 
end do; 
convert(P, list);
end proc:

Example of use:

for j from 1 to 15 do
RandWalk(1500);  L[j]:=nops(%);
print(%%, L[j]);
od:
L:=convert(L, list);  # The list of the lengths of these random walks

Edit.

First you can specify every piece of your space curve by plots[spacecurve]  command then all together by plots[display]  command.

Example:
L1:=plots[spacecurve]([cos(t), sin(t), t], t=0..2*Pi, color=red):
L2:=plots[spacecurve]([cos(t), t-2*Pi, sin(t)+2*Pi], t=2*Pi..4*Pi, color=red):
plots[display](L1, L2, axes=normal, scaling=constrained);

Of course quite easy to write a special procedure that automates this process, if necessary.

Addition.  Here is the procedure which automates the plotting of a piecewise curve in 3D:

PiecewisePlot:=proc(P::piecewise,Opt::list)
local n, m, t, i, L;
uses plots;
n:=nops(P);
if n::odd then error `nops(P) should be even` fi;
t:=op(indets(op(1,P)));
for i from 1 to n/2 do
L[i]:=spacecurve(op(2*i,P),t=op([2*i-1,1,1],P)..op([2*i-1,2,2],P), op(Opt));
od;
display(seq(L[i], i=1..n/2));
end proc:

Example of use:

L:=piecewise(0 <= t and t < 2*Pi,[cos(t), sin(t), t], 2*Pi <= t and t <= 4*Pi, [cos(t), t-2*Pi, sin(t)+2*Pi], t>4*Pi and t<=5*Pi, [1, t-2*Pi, 2*Pi], t>5*Pi and t<=6*Pi, [cos(t),3*Pi,sin(t-5*Pi)+2*Pi]);
PiecewisePlot(L, [color=red, thickness=3, scaling=constrained, axes=normal]);

Edited.
 

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