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

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.