@tomleslie If we use this code for equaion F(x) it is running 

for k from 0 to n do IBVC1 := select(has, CO, b[k]); slv := dsolve({coeff(DE1, p, k), op(IBVC1)}); b[k] := unapply(rhs(slv), x) end do;
F(x) = F(x)+O(p^(n+1));
plot(eval(F(x), p = 1), x = 0 .. 5);
plot(eval(diff(F(x), x), p = 1), x = 0 .. 5);

@Carl Love 

Excellent 

@Kitonum 

Sir in the file attached how we have to put line +symbol in plot 

@Carl Love 

Thanks a lot , excellent sir for quick response

@Carl Love 

Thanks for quick response but still error.

Sir still i am getting an error, may be some problem in  my codes. Please see the codes

restart:
with(plots):
n:=0.75:Eh:=100:mn:=1:t0:=0.2:
a1:=(mn/t0)^((n-1)/(n))*(tb)^(1/n):a3:=Eh/tb:a4:=(Eh)^2:
U1:=(a1/Eh)*(-1+(a3*r/Eh))^(1/n)*(n*a4/(1+n))*(1/a3-(r/Eh))-a1*(-1+(1/tb))^(1/n)*(n/(n+1)*(tb-1)):
plot([seq(seq(eval(U1,tb=j),j in[0.8,0.9,1.0]))],r=0..1,legend = [tb =0.8, tb=0.9,tb =1.0],  labels = ["z ", "U"], labeldirections = ["horizontal", "vertical"],  linestyle = [solid,dash,dot],color = [black, red,green]);
 

@Kitonum 

Thanks a lot

@tomleslie 

thanks a lot

@Carl Love 

still there is an error

restart:
with(plots): 
G:=1300: beta:=0.0075:epsilon:=4.9*10^(-6):Pa:= 100:alpha:=24:Pv:=97:phi:=0.4:

pde1 := G*(beta^2*u(r,z)*diff(u(r,z),z)+v(r,z)*diff(u(r,z),r)) = -diff(p(r,z),z) +beta^2*diff(u(r,z), z$2)+diff(u(r,z),r$2)+(1/r)*diff(u(r,z),r):

pde2 :=G*(beta^2*u(r,z)*diff(v(r,z),z)+v(r,z)*diff(v(r,z),r)) = -diff(p(r,z),r) +beta^2*diff(v(r,z), z$2)+diff(v(r,z),r$2)+(1/r)*diff(v(r,z),r)-(v(r,z)/r^2):

pde3:=beta^2*diff(u(r,z),z)+diff(v(r,z),r)+(v(r,z)/r)=0:

IBCu:={D[1](u)(0, z) = 0,phi*(D[1](u))(1, z)+u(1,z) = 0}:# IBC for u

IBCv:={v(0,z)=0,v(1,z)=epsilon*(p(1,z)+(Pa/alpha)-1)}:# IBC for v

IBCp:={p(r,-1)=0,p(1,z)=(Pv-Pa)/alpha}:# IBC for p

IBC := IBCu union IBCv union IBCp:
PDE:= [pde1,pde2,pde3]:

pdsolve({PDE, IBC}, u(r,z));  

    
num_sol := pdsolve(PDE, { IBC}, numeric);
 

@Carl Love 

As per your suggestion i have make some correction but still error is coming. Even i have seen your code also

https://www.mapleprimes.com/questions/200670-How-Do-I-Plot-Multiple-Pde-Numerical 

restart:
with(plots):
G:=1300: beta:=0.0075:epsilon:=4.9*10^(-6):Pa= 100:alpha:=24:Pv:=97:

pde1 := G*(beta^2*u(r,z)*diff(u(r,z),z)+v(r,z)*diff(u(r,z),r)) = -diff(p(r,z),z) +beta^2*diff(u(r,z), z$2)+diff(u(r,z),r$2)+(1/r)*diff(u(r,z),r):
pde2 :=G*(beta^2*u(r,z)*diff(v(r,z),z)+v(r,z)*diff(v(r,z),r)) = -diff(p(r,z),r) +beta^2*diff(v(r,z), z$2)+diff(v(r,z),r$2)+(1/r)*diff(v(r,z),r)-(v(r,z)/r^2):
pde3:=beta^2*diff(u(r,z),z)+diff(v(r,z),r)+(v(r,z)/r)=0:
PDE:= unapply({pde||(1..3)}, phi):

IBC:=(D[1](u))(0, z) = 0,v(0,z)=0,phi*(D[1](u))(1, z)+u(1,z) = 0,v(1,z)=epsilon*(p(1,z)+(Pa/alpha)-1),p(r,-1)=0,p(1,z)=(Pv-Pa)/alpha:

phiList:= [0.0, 0.15, 0.4]:
Colours:= table(PrList =~ ["blue", "red", "green"]):
for phi in phiList do
     pds:= pdsolve(PDE(phi), IBC, numeric);
     Plots[phi]:= pds:-plot[display](
          u(z,r), r=0.5, linestyle= "solid", colour= Colours[phi],
          legend= sprintf("phi=%2.2f", phi),
          title= "Velocity Profile", labels= ["z", "u"]
     )
end do:
plots:-display([seq(Plots[phi], phi in phiList)]);
 

upload the worksheet

@Christian Wolinski 

I want to plot the curves for different values of m=0.1,0.2,0.3,0.4,0.5(where m is a parameter) 

@Rouben Rostamian  alpha is a parameter. every thing is ok simply how we can plot curves by changing alpha values on 2d plot.

@Rouben Rostamian   thanks 

can you help me to plot the curves for   different values of M = 0.1, 3.0 and 5.0.  in my earlier post where @ tomleslie 5348 helpme to solve the difference scheme . 

maple post link is:crank nicolson system of difference scheme error 

and i am attaching the file .

please help me to plot the curves as mentioned in the post.

bigSys_(2).mw

@Rouben Rostamian  

actually i want to plot the solution in my codes ( where i have used difference scheme) instead of surf data , how to plot those curves .

@Rouben Rostamian  thanks but i want to plot the solution at time t=0.2 for various values of alpha  instead of using surfdata(plotdata,axes=boxed, labels=[x,"u"]): 

in my codes .