In Maple 2015 and 18 I get if I begin with a restart:
restart;
eul := proc(x,y)
    for j from 1 to 3 do
          y[j]:= y + h*fxy(x,y):
          x:= 0 + j*h:
          y := y[j]:
          yt[j] := y:
    end do:
    return yt:
end proc:
##Warnings, but OK
yv := eul(x,y);
#Output yv:=yt
# Then try
y;
The error is "Error, too many levels of recursion"
If you try
eval(yv);
you will find the result:  table([1 = y, 2 = y, 3 = y])
## Now the following will give you the reported error "Error, (in eul) illegal use of a formal parameter" for the obvious reason that you will attempt to assign to 0 the value 0+j*h (in the loop), and that furthermore you will attempt to assign to -2 the value -2[1]; (again in the loop, j=1).

eul(0,-2);
#####################
Notice that Maple actually do use procedures which assign to formal parameters, see the following example in 
?dsolve,numeric,IVP

solproc := proc(N, x, Y, YP)
   YP[1] := Y[2];
   YP[2] := Y[3];
   YP[3] := (2*Y[2] - x*Y[1])/3;
end proc;


 

restart;
f:=x^2*exp(-1.2*x);
res:=maximize(f, x=0..100, location);
a:=subs(indets(res[2],equation),x);

I tried your worksheet and had problems with it. I hate 2d math so converted to 1d math input.
I changed . to * and got for the second derivative w.r.t. theta (called Phi2) with the assignments you have

Phi2 := Phi0*sin((3/2)*Pi*(x-6))*cos(theta)*cos(2*Pi*z)*sin(omega*t);
int(Phi2,z=-1/4..1/4); #Use int: evalf(Int(...)) is for numerical integration which doesn't make any sense when several names don't have numerical values.
Your integral is basically just an integral of cos(2*Pi*z) over z = -1/4..1/4. The rest are just constant factors.

Obviously, u[i]=0 all i is a solution. So I guess the question is if there is another solution.
Anyway, here is a formulation of the problem for n=3.
Example:

P:=LinearAlgebra:-RandomMatrix(3);
U:=Vector(3,symbol=u);
U2:=U^~2;
U=P.U2;
Equate(U,P.U2);
fsolve(%);

As already mentioned by Markiyan this output is due to numerical roundoff.
The exact solution is x(t) = exp(-t), which obviously never becomes zero.

Your equation
ode:=diff(x(t),t)=-abs(x(t));
reduces to diff(x(t),t) = -x(t) on any interval on which x(t) > 0. Since x(0) = 1 > 0, the solution will be x(t) = exp(-t) as long as exp(-t) > 0, which it is for all t.

It is interesting that dsolve (in Maple 2015) returns 2 solutions, but one of those (x(t) = exp(t) ) is obviously wrong!

dsolve({ode,x(0)=1});
combine([%]);
                [x(t) = exp(t), x(t) = exp(-t)]

I checked in Maple 12. It gave only the correct solution.

Leaving out the initial condition in Maple 2015 gives
dsolve(ode);
                      exp(-t)                   
               x(t) = -------, x(t) = exp(t) _C1
                        _C1                     
This is only OK if interpreted correctly: the first _C1 is positive and the second _C1 is negative (or zero). Thus the C's should have had assumptions as are used in outputs from solve.

Doing the same in Maple 12 returns a much better result on implicit form and also better when using the option 'explicit'.
In Maple 2015 one could try using the option 'implicit':
dsolve(ode,implicit); #As in Maple 12
dsolve({ode,x(0)=1},implicit); #Using x(0)=1, but not entirely! The piecewise should have been simplified.

I shall submit an SCR about dsolve({ode,x(0)=1}).

Take a look at
?LinearAlgebra[RowOperation]

if you haven't done so already.

Try this simple example and take a look at the output (besides the warning):

p:=proc(x) y:=x+1 end proc;

#You notice that the output is actually the following, where y is declared local, something you ought to have done yourself. If you actually want y to be a global variable, then replace the word 'local, by 'global'.

p:=proc(x) local y; y:=x+1 end proc;

The system is inconsistent:

restart;
pde1:=diff(u(x, y), x) =-(2/3)*(3*h^3*nu+9*h^2*nu*y-12*nu*y^3+36*x^2*y+56*y^3)/h^3;
pde2:=diff(v(x, y), y) = (2/3)*(36*nu*x^2*y+56*nu*y^3+3*h^3+9*h^2*y-12*y^3)/h^3;
pde3:=diff(u(x, y), y)+diff(v(x, y), x) =-(6*(1+nu))*x*(h^2-4*y^2)/h^3;
pdsolve({pde1,pde2,pde3}); #Notice the warning
#Now by hand:
eq1:=map(int,pde1,x)+(0=f1(y));
eq2:=map(int,pde2,y)+(0=f2(x));
eval(pde3,{eq1,eq2});
diff(%,x,y);
simplify((rhs-lhs)(%)); #Not identically zero

Since you don't give the initial and boundary conditions the following will in part be an example.

restart;
pde:=diff(c(x,t),t)=diff(c(x,t),x,x)-diff(c(x,t),x)-c(x,t);
pdsolve(pde); #Variables separated
pdsolve(pde,build); #Variables separated and the odes solved
pdsolve({pde,c(0,t)=0,c(1,t)=0}); #Example boundary conditions
evalc(%);
about(_Z1); #_Z1 may be some other _Zn if you repeat the code without a restart
#Thus the solution will be of the form
C:=exp(x/2-5*t/4)*Sum(b(n)*sin(n*Pi*x)*exp(-n^2*Pi^2*t),n=1..infinity);
#Now the intial condition c(x,0)=f(x) (whatever it is) must be met:
f(x)=eval(C,t=0);
##
## The solutions give us the idea of making a simple change of dependent variable, which turns 'pde' into the heat equation:
PDEtools:-dchange({c(x,t)=u(x,t)*exp(x/2-5*t/4)},pde,[u]);
op(solve(%,{diff(u(x,t),t)}));
## Thus you can use tricks and methods from textbook examples.

You could use piecewise:

restart;
k := 10; m := 1; g := 10; mu := .2;
Xr := U*t-x(t);
ode1 := m*diff(x(t), t, t) = k*Xr-piecewise(mu*m*g<k*Xr,mu*m*g,0);
U:=1:
res:=dsolve({ode1,x(0)=0,D(x)(0)=0},numeric);
plots:-odeplot(res,[t,x(t)],0..15);
plots:-odeplot(res,[t,diff(x(t),t)],0..15);

The right hand side of ode2 can be written as -mu*m*g*signum(v)-m*A*omega^2*cos(omega*t).
That expression is discontinuous in v at v=0. That becomes a problem when v(t) hits zero (which it does with your initial conditions). When v(t) hits zero the expression changes with the amount 2*mu*m*g.

You could try using dsolve with events:

odeE:=diff(v(t),t)=-0.8*9.8*(2*b(t)-1)-cos(t)^2;
res:=dsolve({odeE,ode1,v(0)=0.25,x(0)=0.5,b(0)=1},numeric,discrete_variables=[b(t)::boolean],events=[[v(t)=0,b(t)=0]]);
plots:-odeplot(res,[t,v(t)],0..1,thickness=3);
plots:-odeplot(res,[t,x(t)],0..1,thickness=3);

The introduction of a floating point number like 1.0 is contagious. Try:
Pi*1.0;
#You will see that you don't get Pi and you don't get its infinite (!) decimal expansion either (obviously).
# sin(Pi*1.0) is going to give you the same as
Pi*1.0;
sin(%);
##To get 0 (well in fact the float 0.) you can do
evalf[11](sin(Pi*1.0));
# assuming that you are using Digits=10 (the default). Any number higher than 11 will do too.

Here is an example using only 3 points. You can just add 3 more points.
Before the animation command I give an image at time t=0:

restart;
with(plots):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]];
ptsL:=[pts[],pts[1]];
pointplot3d(eval(pts,t=0),symbol=solidsphere,symbolsize=50); p1:=%:
pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black); p2:=%:
display(p1,p2);
#Now animating from t = 0 to t = 50:
animate(display,['pointplot3d'(pts,symbol=solidsphere,symbolsize=50),'pointplot3d'(ptsL,connect,thickness=3,color=black)],t=0..50,frames=100);
##Notice the unevaluation quotes around pointplot3d in the animate command. They are there to prevent premature evaluation, i.e. that pointplot3d sees pts and ptsL as containing a symbol t, thus it cannot plot anything.

Assuming that you want real solutions (of which there are infinitely many in all 3 cases) you can use fsolve as follows.
Notice that because cos and cosh are both even, and tan and tanh are both odd, we can restrict ourselves to the positive real axis.

restart;
eq1:=cos(x)*cosh(x)=1;
plot((lhs-rhs)(eq1),x=0..15,-1..1);
fsolve(eq1,x=4..5); #Finds a root between 4 and 5
fsolve(eq1,x=14..15);
eq2:=cos(x)*cosh(x)=-1;
plot((lhs-rhs)(eq2),x=0..15,-1..1);
fsolve(eq2,x=10..12);
eq3:=tan(x)=tanh(x);
plot((lhs-rhs)(eq3),x=0..15,-1..1,discont=true);
fsolve(eq3,x=9..11);
########### solve ########
If you try this version of solve:
res:=solve(eq1,x,AllSolutions);
##
Then you will see that the result is a RootOf-expression:
RootOf(_Z-(2*I)*Pi*_Z5+2*arccosh(1/cos(_Z))*_B4-arccosh(1/cos(_Z)))
You can say that eq1 has just been rewritten some.
You will never get closer to an actual formula for all the roots.
If we use other notation we can write the argument to RootOf as
expr:=x-(2*I)*Pi*n+2*arccosh(1/cos(x))*b-arccosh(1/cos(x));
## _Z5 (or a different number than 5) stands for an arbitrary integer. _B4 can be 0 or 1.
## I have used n and b for those two. Instead of the unknown _Z I have used the name you use in eq1, x.
plot(eval(expr,{n=0,b=0}),x=0..15); #An example
fsolve(eval(expr,{n=0,b=0}),x=4..5); #Compare with the root found earlier

Looking at
indets(C21,name);
after running your worksheet P2.mw I found that theta occurs as just that and also as theta[1] and theta[2].
Furthermore r occurs as just r, but also as r[x], r[y], and r[z].

In general this is a bad idea in Maple because of the special meaning attached to (e.g.) theta[1].

I tried changing theta[1] to theta1, r[x] to rx, and the other 3 similarly.
Not being a friend of 2d-input I copied the whole thing to a new worksheet and converted to 1d input. There I could make the changes with Find/Replace in the Edit menu. (If anybody else is going to try that remember also to replace in the same way semicolon by colon, since the outputs generated are quite large and may force you to shut Maple down).
I didn't receive any errors when running the entire new worksheet.