There are quite a few tools in the package StringTools:

?StringTools

Continuing where I left above saying that the difference in results merits looking into.

There are at least two solutions for n=1 and N=20, and I think for other values as well.

The nice thing about n=1 is that the system is so simple.

Again for completeness I bring the whole code. I don't assign to n or N, but use eval.

restart;
Digits := 15;
with(plots):
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;
sys:=eval({Eqn1, Eqn2, Eqn3},n=1); #Relatively simple
res1:=dsolve(sys union eval({bcs1, bcs2},N=20), numeric); #Using N=20
odeplot(res1,[eta,U(eta)],thickness=3); p1:=%:
odeplot(res1,[eta,V(eta)],thickness=3);
odeplot(res1,[eta,W(eta)],thickness=3);
##Now using an appropriate approximate solution to find the second solution:
res2:=dsolve(sys union eval({bcs1, bcs2},N=20), numeric,
    approxsoln=[U(eta)=eta/20*exp(-eta),V(eta)=-tanh(eta/2),W(eta)=-0.9*tanh(eta/2)]);
odeplot(res2,[eta,U(eta)],thickness=3,color=blue); p2:=%:
odeplot(res2,[eta,V(eta)],thickness=3,color=blue);
odeplot(res2,[eta,W(eta)],thickness=3,color=blue);
display(p1,p2,caption="Two solutions for U");


If this is a problem from some applied area, then the second solution is probably the desired one. I'm guessing that N=20 is a replacement for N=infinity.

Comment. If you use bcs3:=D(U)(N) = 0, D(V)(N) = 0;  instead of bcs2 then the red solution is excluded.

@Dayana It seems that the statement in the help page about initmesh is incorrect:

?dsolve,numeric,bvp

saying (quoting the entire statement):

" 'initmesh'= integer
Integer value that determines the number of points dsolve uses to compute the initial solution profile. In some cases, the default initial 8 point mesh does not have sufficient resolution to obtain the initial solution profile, so increasing this value can give a solution when the default value does not. Its value must be between 8 and 8192."
(My emphasis)

Consider the following simple bvp-problem (whose solution is just y(x) = cos(x) ).
Compare the different usages reported.

restart;
ode:=diff(y(x),x,x)+y(x)=0;
res:=CodeTools:-Usage(dsolve({ode,y(0)=1,D(y)(Pi)=0},numeric,initmesh=30000));
     memory used=86.52MiB, alloc change=83.83MiB, cpu time=3.67s, real time=3.67s, gc time=46.88ms
restart;
ode:=diff(y(x),x,x)+y(x)=0;
res:=CodeTools:-Usage(dsolve({ode,y(0)=1,D(y)(Pi)=0},numeric,initmesh=8192));
     memory used=26.82MiB, alloc change=47.97MiB, cpu time=516.00ms, real time=516.00ms, gc time=0ns
restart;
ode:=diff(y(x),x,x)+y(x)=0;
res:=CodeTools:-Usage(dsolve({ode,y(0)=1,D(y)(Pi)=0},numeric)); #Default
     memory used=4.87MiB, alloc change=32.00MiB, cpu time=47.00ms, real time=46.00ms, gc time=0ns

I think that we must conclude that the statement "Its value must be between 8 and 8192." is incorrect at least all the way back to Maple 12, which is the oldest I have available at my present location.

Also I shall submit an SCR asking for a help page correction (update).
There is a similar statement about maxmesh, but doing the same experiment with maxmesh instead of initmesh doesn't show any significant difference (try several times). That is to be expected since maxmesh is the maximal mesh allowed and if the mesh doesn't need to be raised above the given maxmesh then it isn't.

Thanks for your persistence in your inquiries!

@Dayana When getting the error 'initial Newton iteration not converging' you could try using the option 'initmesh'= N, where N is between 8 and 8192. (8 is the default according to the help for dsolve,numeric,bvp).

Another useful option is to give an approximate solution in the form 'approxsoln'= [U(x)=..., V(x)= ..., W(x)=...].
You can try any expression depending on (in this case) eta.

A third option to consider is the absolute error in the form 'abserr'= eps, where eps is some small positive number. The default is 10^(-6).

The continuation method as you mention is more difficult to use because it involves placing a parameter in a "reasonable" place.

As far as using a midpoint method is concerned, dsolve will often complain itself about problems at an endpoint and suggest using a midpoint method.

Yes, your experience with increasing Digits is one I have had often in bvp-problems. It is a wild goose chase. Generally if I have Digits at 15, I give up raising it because the goose is always going to be ahead of me.
So the solution is generally to make use of other options.

In the present case, I quite agree with I_Mariusz that you should not have right endpoint at 20, but considerably lower. He chose 5, I had success also with 7.
The problems appear because U and V become roughly constant so that U'(eta) = V'(eta) = 0 (almost). When solving for the highest derivatives the quantity U'(eta)^2 +V'(eta)^2 appears in the denominators. That will cause numerical problems although the numerators also will be approaching zero (and I think at a higher rate).

@vv I must have been working on the expanded version of my answer while you posted your reply.
Yes, even for that concrete value of y we don't get a definite answer.

I copied your code into Maple 2015.2 and it ran without any error.
After that I did:
seq(P[j][1,1],j=1..n-1);

and I received as expected 1,1,1.

Did you run the code you showed us just after a restart?
If not try that.

That should be possible in Statistics:-Fit.

?Statistics:-Fit

@roya2001 I don't see any way to get an answer that looks like your approximate solution. That obviously doesn't mean that such a solution doesn't exist, however.

######## An observation:
Let S and G be the expressions for your approxsoln for s(x) and g(x):

S:=cosh(upsilon*x)-cos(upsilon*x)-(cosh(upsilon)+cos(upsilon))*(sinh(upsilon*x)-sin(upsilon*x))/(sinh(upsilon)+sin(upsilon));
G:=sin(((2*n+1)*(1/2))*Pi*x);
##Now one of your boundary conditions is D(g)(1)+1/2*(D(s)(1))^2=0.
#We find
eval(diff(G,x),x=1); # Zero
eval(diff(S,x)^2/2,x=1); #46.0533267839276

Thus that boundary condition is very far from being satisfied by your proposed approximate solution.

Besides what J4James mentions, you also need to enclose the sequence of pdes in {}:
{PDE1, PDE2, PDE3}

@roya2001 The new error is clearly due to the fact that no successes were encountered, thus res is still just a name, not tyhe name of a table.
What did you change since the last version? I cannot keep track of it.

@fbackelj If you are the teacher then it seems to me to be your job to explain the situation to your students after having read the help for LinearAlgebra:-ConditionNumber. It uses the word "estimates" and that surely indicates that you are getting an estimate of the condition number. What are your students using the condition number for? Do they need the exact value for anything?

@Thomas Richard Clearly there could be a bug somewhere, but if I understand the statement in the link given by acer correctly, then there need not be.

The statement in https://www.nag.co.uk/numeric/fl/nagdoc_fl24/html/F07/f07agf.html

says:

"The computed estimate RCOND is never less than the true value ρ, and in practice is nearly always less than 10ρ, although examples can be constructed where RCOND is much larger."

Since the value returned (RCOND) is an estimate of 1/c where c is the (true) condition number, I understand this statement as saying that est_c <= c (always), est_c > c/10 (almost always), but that there are examples where est_c << c/10.

(By est_c I mean the estimated condition number).

####### You may try the following and notice (as Markiyan did) that most often there is agreement, but fairly often est_c is less than c (and not because of rounding).

restart;
eps:=0.1;
N:=0;
R:=Vector():
to 1000 do
use LinearAlgebra in
  A:=RandomMatrix(4,generator=-9.0..9.0);
  Ec:=ConditionNumber(A);
  c:=Norm(A)*Norm(A^(-1))
end use:
  if (c-Ec)/c>eps then
    N:=N+1;
    R(N):=Ec/c;
    print(Ec,c,'ratio'=R(N))
  elif (c-Ec)/c<-10^(-Digits+3) then
    print("Shouldn't happen",Ec,c)
  end if;
end do:


In a trial run I found that in 12 percent of the cases (c-Ec)/c > 0.1 and no "Shouldn't happen".
The ratio min and max
(min,max)(R);
were roughly 0.4 and 0.9.
This is in full compliance with the cited accuracy statement.

After answering and commenting your question at

http://mapleprimes.com/questions/207601-Remove-Assumptions-

I got think that one reason for the statement "This functionality is not intended for end users" may be that complication with assumed variables when saving as an .m file instead as an .mpl file.

Below I was trying to answer your question as formulated.
If you use an .mpl extension instead all your problems are gone. There is nothing to remove.