Dear Maple user

I want to add caption for all four graphs with a,b,c,d at the bottom centre of each graph that should be arrange in 2X2 array and the array bottom caption as Figure 6, also i want to save the output of the array graphs in desktop with jpeg format. sample graph i have uploaded.

restart:
with(plots):
A1:= plot(sin(x), x = -Pi .. Pi,title = "Graph at M^2=1",caption = "Fig a"):# This caption is at the centre bottom of first graph and the title of A1 is M*M (M^2)should be display
A2:= plot(cos(x), x = -Pi .. Pi,caption = "Fig b"):# This caption is at the centre bottom of second graph
A3 := plot3d(1.5^x*sin(y), x = -3 .. 3, y = -3 .. 3,caption = "Fig c"):# This caption is at the centre bottom of Third graph
A4 := plot3d(1.3^x*sin(y), x = -1 .. 2*Pi, y = 0 .. Pi,caption = "Fig d"):# This caption is at the centre bottom of forth graph
display(Array(1 .. 2, 1 .. 2, [[A1, A2], [A3, A4]]), aligncolumns = [1],caption = "Fig 6"):# This caption at the bottom centre of array graphs 

Dear Maple user I am facing error while running the codes  to plot the graph for two data sets .

I am attaching the files.

Error_in_Display_1.mw

Dear maple user i have a partial differential equation whose first approximation solution is w(r,z)=(r^2-h^2)/4 , i want to calculate the two approximations. I am attaching the file please see  

# I dont want to solve directly.

HP_Help.mw

dear maple user help me to correct this error . 

I want to solve the coupled pdes  plot the graphs but unable to get it. please help me. Here are the maple codes

restart:
L:=0.4:
d1:=0.3: d2:=0.3: 

h:=z->piecewise( z<=d1,   1,
                 z<=d1+L,   1-(delta/(2))*(1 + cos((2*(Pi)/L)*(z - d1 - L/2))), 
                 z<=d1+L+d2,  1 ):

a:=x4*S*Gr/(4*x1*x5)*sin(alpha):
b:=1/Da+(x3*M/(x1*(1+m^2))):
c:=(1/x1)*DP:
x1:=1/((1-phi1)^2.5*(1-phi2)^2.5):
x5:=khnf/kf:
x3:=shnf/sf:
x4:=(1-phi2)*((1-phi1)+phi1*(RBs1)/(RBf))+phi2*RBs2/RBf:
shnf:=sbf*((ss2+2*sbf-2*phi2*(sbf-ss2))/(ss2+2*sbf+phi2*(sbf-ss2))):
sbf:=sf*((ss1+2*sf-2*phi1*(sf-ss1))/(ss1+2*sf+phi1*(sf-ss1))):
ss2:=2.7*10^(-8):
ss1:=59.6*10^(6):
sf:=6.67*10^(-1):
khnf:=kbf*((ks2+2*kbf-2*phi2*(kbf-ks2))/(ks2+2*kbf+phi2*(kbf-ks2))):
kbf:=kf*((ks1+2*kf-2*phi1*(kf-ks1))/(ks1+2*kf+phi1*(kf-ks1))):
ks1:=401:
ks2:=76.5:
kf:=0.4972:
RBs1:=8933*16.7*10^(6):
RBs2:=6320*18*10^(6):
RBf:=1063*1.8*10^(6):
DP:=(1/192)*(18432*F-77*b^2*h(z)^8+308*b*h(z)^8*a+768*h(z)^6*a+12*b*h(z)^8)/(h(z)^4*(b*h(z)^2-6)):
conds:= w(h(z),0)=0, theta(h(z),0)=0, D[1](w)(r,0)=0, D[1](theta)(r,0)=0:
pdes:= diff(w(r, z), r, r)+(1/r)*(diff(w(r, z), r))-(1/Da+(x3/x1)*M)*w(r, z)-(1/x1)*(DP+theta*x4*sin(alpha))=0,
         diff(theta(r, z), r,r)+(1/r)*(diff(theta(r, z), r))+(S/x5)=0:
pars:= { F=1.5,Da=0.1,phi1=0.01,phi2=0.02,alpha=Pi/4,S=0.1,Gr=2,m=0.5}:
         
  PrVals:=[2, 5, 7,9]:
  colors:=[red, green, blue, black]:
  for j from 1 to numelems(PrVals) do
      pars1:=`union`( pars, {M=PrVals[j]}):
      pdSol:= pdsolve( eval([pdes], pars1),
                       eval([conds], pars1),
                       numeric
                     );
      plt[j]:=pdSol:-plot( w(r,z), z=1, r=0..1, numpoints=200, color=colors[j]);
  od:
  plots:-display( [seq(plt[j], j=1..numelems(PrVals))]);
Error, final value in for loop must be numeric or character
Error, unable to execute seq
 PrVals:=[2, 5, 7,9]:
  colors:=[red, green, blue, black]:
  for j from 1 to numelems(PrVals) do
      pars1:=`union`( pars, {M=PrVals[j]}):
      pdSol:= pdsolve( eval([pdes], pars1),
                       eval([conds], pars1),
                       numeric
                     );
      plt[j]:=pdSol:-plot( theta(r,z), z=1, r=0..1, numpoints=200, color=colors[j]);
  od:
  plots:-display( [seq(plt[j], j=1..numelems(PrVals))]);
Error, final value in for loop must be numeric or character
Error, unable to execute seq

 PrVals:=[2, 5, 7,9]:
  colors:=[red, green, blue, black]:
  for j from 1 to numelems(PrVals) do
      pars1:=`union`( pars, {M=PrVals[j]}):
      pdSol:= pdsolve( eval([pdes], pars1),
                       eval([conds], pars1),
                       numeric
                     );
      plt[j]:=pdSol:-plot( DP,  z=0..1, numpoints=200, color=colors[j]);
  od:
  plots:-display( [seq(plt[j], j=1..numelems(PrVals))]);

Dear Maple users Help me to  get the desire graph for this codes. 

restart:
with(plots):
with(IntegrationTools):
h:=z->piecewise( z<=d+1,   1,
                 z<=d+4,   1-(delta/(2))(1 + cos(2(Pi)*(z - 1 - 1/2))),                                                           z<=d+6,   1 ):
w0:=(-c*h(z)^2/4)+(3/64)(b*c-4*a)*h(z)^4+(19/2304)*b(b-4*a)*h(z)^6:
w1:=(c/4)+(1/16)*(4*a-b*c)*h(z)^2:
w2:=(1/256)(4(b*c-4*a)-b*h(z)^2):
w3:=(1/2304)b(b-4*a):
a:=(x4*S*Gr)*sin(alpha)/(4*x1*x5):
b:=(1/Da)+(x3*M/(x1*(1+m^2))):
c:=(1/x1)*Dp:
Dp:=96*x1/((6-b*h(z)^2)h(z)^4)(F+(a*h(z)/24)-((11/6144)b(b-4*a)*h(z)^8)):
x1:=1/((1-phi1)^2.5*(1-phi2)^2.5):
x2:=(1-phi2)((1-phi1)+phi1*Rs1/Rf)+phi2(Rs2/Rf):
x3:=(shnf)/(sf):
x4:=(1-phi2)((1-phi1)+phi1(RBs1)/(RBf))+phi2*((RBs2)/(RBf)):
x5:=khnf/kf:
shnf:=sbf*((ss2+2*sbf-2*phi2*(sbf-ss2))/(ss2+2*sbf+phi2*(sbf-ss2))):
sbf:=sf*((ss1+2*sf-2*phi1*(sf-ss1))/(ss1+2*sf+phi1*(sf-ss1))):
khnf:=kbf*((ks2+2*kbf-2*phi2*(kbf-ks2))/(ks2+2*kbf+phi2*(kbf-ks2))):
kbf:=kf*((ks1+2*kf-2*phi1*(kf-ks1))/(ks1+2*kf+phi1*(kf-ks1))):
RBs1:=(8933*16.7*10^6):
RBf:=(1063*1.8*10^6):
RBs2:=6320*18*10^6:
kf:=0.492:
sbf:=6.67*10^(-1): ss2:=2.7*10^(-8):
sf:=6.67*10^(-1):ss1:=59.6*10^(6):
ks2:=76.5:kf:=0.492: ks1:=401:
phi1:=0.01: phi2:=0.02:alpha:=Pi/4:m:=0.5:Da:=0.1:Gr:=5:delta:=1:S:=0.5:  d:=1:
                                      
W1:=w0+w1*r^2+w2*r^4+w3*r^6:

by varing M =2,5,7 and r varies from 0 to 1 i want this type of graphs.  please see the sample graphs