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.

In the first session l1 is renamed to the escaped local l1~. Then a=sqrt(l1~) is saved.

In the second session what is read in is a as sqrt(l1~) (actually (l1~)^(1/2)).

You can access that l1~ in (at least) two ways:

L:=op(1,a); #You could use l1 instead of L
subs(L=2,a); #Works

##  alternatively (and this time defining l1 as the old globalized escaped local):

l1:=cat(l1,`~`);
subs(l1=4,a);
simplify(%);

I don't see any evidence that assumptions are stored, though.

## OK one more way without doing any of the other is redefining a:
a:=subsop(1=l1,a);

##################################
Since you may have several assumed variables, you may wish a programmatic approach.
Example:

restart;
assume(l1>0,l2>0,l3>0); #Any assumptions
expr := sqrt(l1)+sin(x*l2)+exp(-l3)*l1;
save expr, "E:/test3.m";
#############
restart;
read "E:/test3.m";
expr;
nms:=convert(indets(expr,name),list);
nmsS:=convert~(nms,string);
L:=StringTools:-Substitute~(nmsS,"~",""); #Removing "~"
L1:=parse~(L);
S:=nms=~L1;
Expr:=subs(S,expr);
eval(Expr,{l1=4,l2=5,l3=6});
simplify(%);

In the help page for LinearAlgebra:-ConditionNumber we find the statement:

For floating-point Matrices, the norm of the inverse is computed by using estimates of the norms of the inverses of the triangular factors of the Matrix obtained by LU or Cholesky factorization.

Thus the method used is different from the one used in the non float case. After a few comparisons this difference doesn't seem to have anything to do with roundoff, but is rather a somewhat different measure of the condition number. Whether this should be considered a bug or just an irrelevant difference I shall leave to others to decide. (After all you might as well have used another norm).

In the exact (non-float) case this is used:

use LinearAlgebra in
Norm(A, infinity)*Norm(MatrixInverse(A), infinity)
end use;

If you repeat that with B you get 6.

If you try

lprint(u);

you find

21.1896000000000*z(t)-.784800000000000*sin(theta(t))*(0.5000000000e-1*2^(1/2)*cos(`thet1(t)`)+0.5000000000e-1*2^(1/2)*sin(`thet1(t)`))+.784800000000000*cos(theta(t))*sin(phi(t))*(-0.5000000000e-1*2^(1/2)*cos(`thet1(t)`)+0.5000000000e-1*2^(1/2)*sin(`thet1(t)`))-.784800000000000*sin(theta(t))*((1/2)*2^(1/2)*(.10*cos(`thet1(t)`)*cos(`thet2(t)`)-.10*sin(`thet1(t)`)*sin(`thet2(t)`)+.2*cos(`thet1(t)`))+(1/2)*2^(1/2)*(.10*sin(`thet1(t)`)*cos(`thet2(t)`)+.10*cos(`thet1(t)`)*sin(`thet2(t)`)+.2*sin(`thet1(t)`)))+.784800000000000*cos(theta(t))*sin(phi(t))*(-(1/2)*2^(1/2)*(.10*cos(`thet1(t)`)*cos(`thet2(t)`)-.10*sin(`thet1(t)`)*sin(`thet2(t)`)+.2*cos(`thet1(t)`))+(1/2)*2^(1/2)*(.10*sin(`thet1(t)`)*cos(`thet2(t)`)+.10*cos(`thet1(t)`)*sin(`thet2(t)`)+.2*sin(`thet1(t)`)))

You will notice that thet1(t) doen't appear as such, but as `thet1(t)`. This means that you have made it into a name probably unintentionally. In the present context, however, this is not so bad an idea. Then you don't even need Physics:-diff, you can use ordinary diff:

diff(u,`thet1(t)`);

You have done the same with thet2(t).

A simple way of changing the original u to the intended version:
w:=evalindets(u,name,parse);
lprint(%);
diff(w,thet1(t)); #Now this gives an error
Physics:-diff(w,thet1(t)); #This gives the result

I used piecewise in this. `dsolve/bvp/convertsys` doesn't seem to like that: It issues an error about not being able to convert to an explicit first order system. See note at the bottom.
After that I tried an ivp problem. dsolve had no problem with that.

restart;
sys:={diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))-piecewise(psi[S](t)>-c,upsilon,0),
diff(psi[S](t), t) = eta*K(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2), diff(psi[Iota](t), t) = eta*S(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2)};

bcs:={S(0)=1,K(0)=1,psi[S](T)=1,psi[Iota](T)=1};

params:={N=1,w=1,eta=1,upsilon=1,c=0.5,T=1};
res:=dsolve(eval(sys union bcs,params),numeric); #Error mentioned above
##The system looks good to me:
eval(sys union bcs,params);
# So try an ivp:
ics:={S(0)=1,K(0)=1,psi[S](0)=1,psi[Iota](0)=5};
res:=dsolve(eval(sys union ics,params),numeric); #No problem
plots:-odeplot(res,[t,psi[S](t)],0..2);
plots:-odeplot(res,[t,S(t)+K(t)],0..2.1,thickness=3);
plots:-odeplot(res,[t,S(t)],0..2);

NOTE: The problem `dsolve/bvp/convertsys` had was apparently caused by the use of S in psi[S]. S is thus appearing in two roles, as the S in S(t) and as the S in psi[S]. Changing the latter to e.g. s (lower case s) makes dsolve work on the bvp.