Don't ever try to execute a whole block of statements without having executed them one after the other and observed the output.
I split your whole block of code into its single commands and started executing from the top. I didn't get very far. After having executed Digits:=13, executing the next one:

sys_ode := diff(ph(t), t) = urj*prj(t)+urd*prd(t)+ugd*pgd(t)-yc*ph(t), diff(pc(t), t) = yc*ph(t)-(yrj+yrd+ygd+yf)*pc(t), diff(prj(t), t) = yrj*pc(t)-urj*prj(t), diff(prd(t), t) = -urd*prd(t)+yrd*pc(t), diff(pgd(t), t) = -ugd*pgd(t)+ygd*pc(t), diff(pf(t), t) = yf*pc(t)*simplify(add(u, u = sys_ode));

resulted in the error message
   
Error, recursive assignment

That error is issued since you have the as yet unassigned sys_ode on the left of := but also on the right buried in add(u, u = sys_ode).

This assignment would make sense, if sys_ode had previously been defined to be something, but it hasn't.
A very simple example of the same thing:

restart;
a:=a+1;
Error, recursive assignment

Do you mean that by R = 2*b*(b^2-9*b-6)/(10*b^3+21*b^2+16*b+4) is defined a function of one variable (b)?

If so, you will notice from the plot below (and in other ways) that R is negative for b=1..5 and not particularly large in absolute value.

restart;
R := 2*b*(b^2-9*b-6)/(10*b^3+21*b^2+16*b+4);
plot(R,b=1..5);
solve(R=0,b);
evalf(%);

restart;
E:=-2*P1*A[Q1]+A[0]^2+A[Q1]^2-2*A[0]*rho*P2/(rho+Q1)+A[Q3]^2-2*P3*A[Q3]+P1^2+rho^2*P2^2/(rho+Q1)^2+P3^2+m^2*c^4;
mtaylor(E,[A[Q1]=P1,A[0]=P2*rho/(rho+Q1),A[Q3]=P3],6);
## You can obtain the same result with:
Student:-Precalculus:-CompleteSquare(E,[A[0],A[Q1],A[Q3]]);

In your second example (a first order system of 8 odes) you have a function, which basically switches from omega[sw] to omega[st] in a very short time. The way you have implemented it, the switch is continuous. I'm assuming that this switch is implemented in an ad hoc manner and could as well be done differently.

Thus I shall present a way to do it discontinuosly by using the events option in dsolve/numeric:

restart;
Digits:=15:
K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);
for i to 4 do
  r:=sqrt((u[i](t))^2+(v[i](t))^2);
  Equ[i]:=diff(u[i](t),t)=Au*(1-r^2)*u[i](t)-omega[i](t)*v[i](t);
  Eqv[i]:=diff(v[i](t),t)=Av*(1-r^2)*v[i](t)+omega[i](t)*u[i](t)+(K. <seq(v[i](t),i=1..4)>)[i];
  EqSys[i]:=[Equ[i],Eqv[i]]
end do:
params:={omega[st]=2*Pi,omega[sw]=4*2*Pi,Au=5,Av=5};

sys:={seq(EqSys[i][],i=1..4)};
##In ic the boolean variable b[i](t) is given initial value 0 if v[i](0) < 0 and 1 otherwise:
ic:={u[1](0)=0, v[1](0)=-.01,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1,b[1](0)=0,b[2](0)=0,b[3](0)=1,b[4](0)=1};
## The b[i]'s take on values 0 and 1 only:
sysE:=subs(seq(omega[i](t)=b[i](t)*omega[st]+(1-b[i](t))*omega[sw],i=1..4),sys);
resE:=dsolve(eval(sysE union ic,params), numeric,maxfun=0,discrete_variables=[seq(b[i](t)::boolean,i=1..4)],
        events=[seq([side(v[i](t),post),b[i](t)=1-b[i](t)],i=1..4)]);
plots:-odeplot(resE,[u[1](t),v[1](t)],0..15,numpoints=2000);
plots:-odeplot(resE,[u[1](t),v[1](t)],100-.6..100,numpoints=2000);
plots:-odeplot(resE,[[t,u[1](t)],[t,v[1](t)]],0..15,numpoints=2000,size=[1800,200]);
plots:-odeplot(resE,[t,u[1](t)],0..15,numpoints=2000,size=[1800,200]);
plots:-odeplot(resE,[[t,b[1](t)],[t,v[1](t)]],0..15,numpoints=2000,size=[1800,200],thickness=[3,1]);
plots:-odeplot(resE,[u[1](t),v[1](t),u[2](t)],0..15,numpoints=2000);
plots:-odeplot(resE,[u[1](t),v[1](t),u[3](t)],10..15,numpoints=2000);
##For an explanation of the trigger side(y(t),post), see:
?dsolve,numeric,events

This works:

odeplot(p,[r(theta)*cos(theta),r(theta)*sin(theta)],labels=[x,y]);

Comments:
odeplot doesn't complain here, but the result is wrong:
odeplot(p,[theta,r(theta)],0..2*Pi,coords=polar,axiscoordinates=polar);

plot works fine with that:

plot(sin(2*theta),theta=0..2*Pi, coords = polar, axiscoordinates = polar);
plot(sin(2*theta),theta=0..2*Pi, coords = polar);

You don't really have an equation (no equality sign), but I shall interpret this in the standard Maple way that you intend f = 0.

restart;
f := y(x)+5*diff(y(x), x)+2*x^3*D(y)(0)+3*y(0)=0;
## Introduce names y0 and y1 instead of y(0) and D(y)(0) (=eval(diff(y(x),x) ,x=0) ):
##
ode:=y(x)+5*diff(y(x), x)+2*x^3*y1+3*y0=0;
res:=dsolve({ode,y(0)=y0});
diff(rhs(res),x);
eq:=eval(%,x=0) = y1;
solve(eq,{y1});
sol:=subs(%,res);
## Answer:
                      y(x) = (8/5)*x^3*y0-24*x^2*y0+240*y0*x-1203*y0+1204*exp(-(1/5)*x)*y0
where y0 can be chosen as you wish.

##Test of answer:
eval(f,y=unapply(rhs(sol),x));
##Outout                0=0

##You may want to factor:
factor(sol);

A:=Matrix([[1,2],[3,4]]);
##Either
A.A^%T;
## or
A.LinearAlgebra:-Transpose(A);

arctan(x)+arctan(x/2)=Pi/2;
solve(%,x);
simplify(%);

p:=3.5*x^2+3.2*x-6.5+88.3*x*y-y^3+a*y;
evalindets(p,float,round);

Using as an approximate solution the result from a previous value of n can help.
I repeat the whole code and I'm using eta=N=20 as the right end point.

restart;
Digits := 15;
with(plots):
#n:=13/10: # n unassigned as of now
mu(eta):=(diff(U(eta),eta)^(2)+diff(V(eta),eta)^(2))^((n-1)/(2)):
Eqn1 := 2*U(eta)+(1-n)*eta*(diff(U(eta), eta))/(n+1)+diff(W(eta), eta) = 0;
Eqn2 := U(eta)^2-(V(eta)+1)^2+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(U(eta), eta))-mu(eta)*(diff(U(eta), eta, eta))-(diff(U(eta), eta))*(diff(mu(eta), eta)) = 0;
Eqn3 := 2*U(eta)*(V(eta)+1)+(W(eta)+(1-n)*eta*U(eta)/(n+1))*(diff(V(eta), eta))-mu(eta)*(diff(V(eta), eta, eta))-(diff(V(eta), eta))*(diff(mu(eta), eta)) = 0;
bcs1 := U(0) = 0, V(0) = 0, W(0) = 0;
bcs2 := U(N) = 0, V(N) = -1;
N:=20:
################
AP:=NULL;
for n from 1 by 0.1 to 1.9 do
  res[n]:= dsolve({Eqn1, Eqn2, Eqn3, bcs1, bcs2}, initmesh = 10000, output = listprocedure, numeric,AP,abserr=1e-5);
  AP:=approxsoln=res[n]
end do;
#error at n=1.9
indices(res,nolist);
##Now going in increments of 0.01 from 1.81 and up to n=1.9.
AP:=approxsoln=res[1.8]: #Known from above
for n from 1.81 by 0.01 to 1.9 do
  res[n]:= dsolve({Eqn1, Eqn2, Eqn3, bcs1, bcs2}, initmesh = 10000, output = listprocedure, numeric,AP,abserr=1e-5);
  AP:=approxsoln=res[n]
end do;
#Apparent success all the way.
indx:=sort([indices(res,nolist)]);
#Animations:
display(seq(odeplot(res[indx[i]],[eta,U(eta)],caption=typeset("n = ",indx[i])),i=1..nops(indx)),insequence);
display(seq(odeplot(res[indx[i]],[eta,V(eta)],caption=typeset("n = ",indx[i])),i=1..nops(indx)),insequence);
display(seq(odeplot(res[indx[i]],[eta,W(eta)],caption=typeset("n = ",indx[i])),i=1..nops(indx)),insequence);
#One of the animations:

The difference of the results for N=5 and N=20 is striking and merits looking into.

You have an equality sign (=) instead of the assignment sign (:=) in Pr=1.

I suggest not using colon (:) but semicolon (;) so that you would see that Pr "survived" in the equations as just Pr.

I had success with initmesh = 1024.

The reason for the error message is that with an extra boundary condition an unassigned parameter (such as Pr) can be determined.

But surely you just made a simple mistake.

In Maple arcsin is the principal branch, i.e. it is the inverse of the sine function restricted to the interval [-Pi/2, Pi/2].

To get what you want you could use solve:

solve({sin(x)=-1,x>=Pi/2,x<=3*Pi/2},x,allsolutions);

which returns  {x = (3/2)*Pi}.

Unfortunately, the more general version

solve({sin(x)=y,x>=Pi/2,x<=3*Pi/2},x,allsolutions);

doesn't return a definite answer. It returns

             {x = -2*arcsin(y)*_B2+Pi*_B2+2*Pi*_Z2+arcsin(y)}

which really is the formula for all solutions, since _B2 stands for 0 or 1 and _Z2 for any integer.
By taking _B2 = 0 and _Z2 = 1 we get that value of x which lies in the interval [Pi/2, 3*Pi/2].
That we have to do that "by hand" is regrettable.
Inadequacies in this area came up quite recently in a question by Kitonum:

http://www.mapleprimes.com/questions/207525-Strange-Bugs-Of-Solve-Command

Let X and Y be the vectors of x and y data and T the vector of time data. Further let xe and ye be the functions giving the exact x- and y-values according to the theoretical model.

Then you can use Fit with stacked data as follows:

t2:=max(T)+1;
T2:=Vector([T,T+~t2]);
XY:=Vector([X,Y]);
##The model in piecewise form:
xyt:=piecewise(t<t2,xe(t),ye(t-t2));
sol:=Statistics:-Fit(xyt,T2,XY,t);
solx:=simplify(sol) assuming t<t2;
soly:=simplify(sol) assuming t>=t2;
soly:=eval(soly,t=t+t2);
plot([solx,soly,t=0..1]);

#This seems to work just fine, but is of course somewhat of a trick.

If X and Y are the vectors of x and y data, and T is the vector of time data, then you could minimize the sum of squares using Optimization:-LSSolve.
The list of residuals would be

[seq(X[i]-x(T[i]),i=1..n),seq(Y[i]-y(T[i]),i=1..n)]

where n = numelems(X) = numelems(Y) = numelems(T) and where x(t), y(t) gives the theoretical values in terms of t and your parameter(s).

Note.

I just tried it on the model you described. I used made up data. It worked.

Another note:

I also tried a nonlinear model using dsolve/numeric with option parameters. There I had success with using NLPSolve on the sum of squares directly (but in fact also with LSSolve).

I can post the code if you are interested.

By doing B:=A;  when A is a Matrix (rtable) you are just making B refer to the same mutable object as A. Thus if one of them is changed, the other is too, because what is changed is what they both refer to.

To create a new Matrix B (as you intended, I presume) use either copy or LinearAlgebra:-Copy.

B:=copy(A);

Read the help pages for copy and LinearAlgebra:-Copy:

?copy
?LinearAlgebra:-Copy

Comment. The important word above is mutable.
If you do:
a:=2;
b:=a;
b:=56;
then surely a still evaluates to 2. What a refers to, is the integer 2, which is not a mutable object (it cannot be changed).
By doing b:=a; you are making b refer to whatever a evaluates to at the time of execution of that command. Thus here b will refer to 2. After that the command b:=56; will, when executed, make b refer to 56 instead of 2.