@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 .

@Rouben Rostamian  

with out adding a parameter  we can plot the graphs. the codes are mentioned below(codes refered from Ralph book) . But if we add a parameter "alpha=0.1,0.5,1,2,5" by taking t=0.2(suppose) how we can plot the graph for different values of alpha .0.1,0.5,1,2,5, the codes are mentioned in post obove

restart:
with(linalg):with(plots): 
ge[1]:=diff(u[1](x,t),t)=diff((u[2](x,t)-1)*diff(u[1](x,t),x),x)+(16*x*t-2*t-16*(u[2](x,t)-1))*(u[1](x,t)-1)+10*x*exp(-4*x):
ge[2]:=diff(u[2](x,t),t)=diff(u[2](x,t),x$2)+diff(u[1](x,t),x)+4*(u[1](x,t)-1)+x^2-2*t-10*t*exp(-4*x): 
bc1[1]:=u[1](x,t)-1: 
bc1[2]:=u[2](x,t)-1: 
bc2[1]:=3*u[1](x,t)+diff(u[1](x,t),x)-3:
bc2[2]:=5*diff(u[2](x,t),x)-evalf(exp(4))*(u[1](x,t)-1):
IC[1]:=u[1](x,0)=1: 
IC[2]:=u[2](x,0)=1: 
NN:=2: 
N:=2:
L:=1:
for i to NN do  
dydxf[i]:=1/2*(-u[2,i](t)-3*u[0,i](t)+4*u[1,i](t))/h: 
dydxb[i]:=1/2*(u[N-1,i](t)+3*u[N+1,i](t)-4*u[N,i](t))/h:
dydx[i]:=1/2/h*(u[m+1,i](t)-u[m-1,i](t)); 
d2ydx2[i]:=1/h^2*(u[m-1,i](t)-2*u[m,i](t)+u[m+1,i](t)):od:
 for i to NN do bc1[i]:=subs(diff(u[1](x,t),x)=dydxf[1],
diff(u[2](x,t),x)=dydxf[2],u[1](x,t) 
=u[0,1](t),u[2](x,t)=u[0,2](t),x=0,bc1[i]):od: 
for i to NN do bc2[i]:=subs(diff(u[1](x,t),x)=dydxb[1],
diff(u[2](x,t),x)=dydxb[2],u[1](x,t) 
=u[N+1,1](t),u[2](x,t)=u[N+1,2](t),x=L,bc2[i]):od:
for i to NN do eq[0,i]:=bc1[i];eq[N+1,i]:=bc2[i]:od: 
for i from 1 to N do eq[i,1]:=diff(u[i,1](t),t)= subs(diff(u[1](x,t),x$2) =
subs(m=i,d2ydx2[1]), 
diff(u[2](x,t),x$2) = subs(m=i,d2ydx2[2]),diff(u[1](x,t),x) =
subs(m=i,dydx[1]),diff(u[2](x,t),x) = subs(m=i,dydx[2]),u[1](x,t)=u[i,1](t), 
u[2](x,t)=u[i,2](t),x=i*h,rhs(ge[1])):od:

for i from 1 to N do eq[i,2]:=diff(u[i,2](t),t)= subs(diff(u[1](x,t),x$2) =
subs(m=i,d2ydx2[1]), 
diff(u[2](x,t),x$2) = subs(m=i,d2ydx2[2]),diff(u[1](x,t),x) =
subs(m=i,dydx[1]),diff(u[2](x,t),x) = subs(m=i,dydx[2]),u[1](x,t)=u[i,1](t),
u[2](x,t)=u[i,2](t),x=i*h,rhs(ge[2])):od: 

for i to NN do u[0,i](t):=(solve(eq[0,i],u[0,i](t))):od:

 for i to NN do u[N+1,i](t):=(solve(eq[N+1,i],u[N+1,i](t))):od:

 h:=L/(N+1): 

for i from 1 to N do eq[i,1]:=eval(eq[i,1]):od: 
 for i from 1 to N do eq[i,2]:=eval(eq[i,2]):od:

eqs:=seq(seq((eq[i,j]),i=1..N),j=1..NN): 
Y:=seq(seq(u[i,j](t),i=1..N),j=1..NN): 

 ICs:=seq(u[i,1](0)=rhs(IC[1]),i=1..N),seq(u[i,2](0)=rhs(IC[2]),i=1..N): 

sol:=dsolve({eqs,ICs},{Y},type=numeric,stiff=true,maxfun=1000000,abserr=1e-6,relerr=1e-5,output=listprocedure):

for j to NN do for i to N do U[i,j]:=subs(sol,u[i,j](t)):od:od: 

for i to NN do U[0,i]:=subs(u[1,1](t)=U[1,1],u[1,2](t)=U[1,2],
u[2,1](t)=U[2,1],u[2,2](t)=U[2,2],u[0,i](t)):od:
 for i to NN do U[N+1,i]:=eval(subs(u[N,1](t)=U[N,1],u[N,2](t)=U[N,2],
u[N-1,1](t)=U[N-1,1],u[N-1,2](t)=U[N-1,2],u[N+1,i](t))):od:
tf:=1.: 
M:=30: 
T1:=[seq(tf*i/M,i=0..M)]: 
PP:=matrix(N+2,M+1): 
for i from 1 to N+2 do PP[i,1]:=evalf(subs(x=(i-1)*h,rhs(IC[1]))):od: 
for i from 1 to N+2 do for j from 2 to M+1 do
PP[i,j]:=evalf(subs(t=T1[j],U[i-1,1](t))):od:od:
plotdata:=[seq([seq([(i-1)*h,T1[j],PP[i,j]],i=1..N+2)],j=1..M+1)]:
t:=0.2:
surfdata(plotdata,axes=boxed, labels=[x,"u"]):
 

@tomleslie 

I want to plot the graphs for U, T. C by changing different parameters 

@tomleslie  thanks a lot for your help. 

i want to plot the graphs for different values  of parameter M=2,4  and fixing t=j=0.2 and   y=i=0 to 4on x axis and U on y axis. I am thankful to you for rectifying the error and help me to plot the graph .

sample plot is attached

@tomleslie 

Thanks for your annotations because of that i have try to  learn how to write the codes.

Please help me  to plot the curves .

I want to plot the 2 dimension plot for different values of M . I have three coupled difference scheme .> i am attaching the codes  

restart; 
# Parameter values:
 Pr:=0.71:E:=1:A:=0:Sc:=0.02: K:=1:

a := 0: b := 1: N := 9:
h := (b-a)/(N+1): k := (b-a)/(N+1):

 lambda:= 1/h^2:  lambda1:= 1/k^2:
# Initial conditions
for i from 0 to N do 
  U[i, 0] := h*i+1:
end do:

for i from 0 to N do 
  T[i, 0] := h*i+1:
end do:
for i from 0 to N do 
  C[i, 0] := h*i+1:
end do:

# Boundary conditions
for j from 0 to N+1 do 
  U[0, j] := exp(A*j*lambda); 
  U[N+1, j] := 0 
end do:

for j from 0 to N+1 do 
  T[0, j] := j*lambda1; 
  T[N+1, j] := 0 
end do:

for j from 0 to N+1 do 
  C[0, j] := j*lambda1; 
  C[N+1, j] := 0 
end do:

# Boundary conditions
#Discretization Scheme
for i to N do 
  for j from 0 to N do 
    eq1[i, j]:= lambda1*(U[i, j+1]-U[i, j]) = (Gr/2)*(T[i, j+1]+T[i,j])+(Gr/2)*(C[i, j+1]+C[i,j])+(lambda^2/2)*(U[i-1,j+1]-2*U[i,j+1]+U[i+1,j+1]+U[i-1,j]-2*U[i,j]+U[i+1,j])-(M/2)*(U[i,j+1]+U[i,j]) 
  end do
end do:

for i to N do 
  for j from 0 to N do 
    eq2[i, j]:= lambda1*(T[i, j+1]-T[i, j]) = (1/Pr)*(lambda^2/2)*(T[i,j+1]-2*T[i,j+1]+T[i+1,j+1]+T[i-1,j]-2*T[i,j]+T[i+1,j])+(E*lambda^2)*((U[i+1,j]-U[i,j])^2) 
  end do
end do:

for i to N do 
  for j from 0 to N do 
    eq3[i, j]:= lambda1*(C[i, j+1]-C[i, j]) = (1/Sc)*(lambda^2/2)*(C[i,j+1]-2*C[i,j+1]+C[i+1,j+1]+C[i-1,j]-2*C[i,j]+C[i+1,j])+(K/2)*((C[i,j+1]+C[i,j])) 
  end do
end do:
#for  different values of M=2,4,6 i want to plot the graph as mentioned above 
 

For the coupled  differential  equations  if we apply GalerkinFEM we can get the solution. Uptill  now I did not see galerkin fem in maple.

@vv 

Because once the above TECHNIQUE is applied to simple problem,  we can applied  to difficulty  problem  also 

@tomleslie 

We can see this question also

https://www.mapleprimes.com/questions/225117-How-To-Solve-This-Error-In-Maple

@tomleslie 

Thanks, actually I want to apply the procedure given in the reference paper , in that paper they use MATLAB to solve  by galerkin technique and written the difference scheme  and plot the graphs.

IS IT POSSIBLE IN MAPLE ALSO CAN WE SOLVE BY THE SAME TECHNIQUE ??????.if it so  then number of non linear coupled equation  can be solve.

Requesting a maple user's to apply the same TECHNIQUE  given in the reference paper.  I have also tried it.

Dear maple user, any one know why it happen, what's wrong with the codes

@acer Thanks a lot for  patiently answer my queries.

Sorry I don't want to  screw your  edits .

Just want to learn maple programming.

One small suggestions i want to know from your side

I want to plot the graph for different values of k=0.1,0.2,0.3,04.

Is it possible to change the last line  of you code to get 4 graphs  in 2X2  matrix form

plot([seq(eval(igRe,gamma2=j),j=[0,0.02,0.06])],gamma1=0.02..0.1,  adaptive=false,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],  linestyle = [solid,dash,dot],
     color = [black,black,black],    labels=[gamma1,'Re(Dp)'],     gridlines=false, axes=boxed);

changed as   

display(Matrix(2,2,[seq(plot([seq(eval(igRe(k),gamma2=j),j=[0,0.02,0.06])],gamma1=0.02..0.1,adaptive=false,legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],linestyle = [solid,dash,dot], color = [black,black,black],
     labels=[gamma1,'Re(Dp)'],gridlines=false, axes=boxed),k in [0.1,0.2,0.3,0.4])]));
while changing  i am getting an error 

@acer  thanks

I am still getting stuck in the plot if i add further terms.  i have add some more terms but while running the codes it will take more times and the system stuck. See what is  wrong in the codes that i  am unable to plot the graphs 

Help1.mw

@acer 

thanks a lot  for making imaginary part as zero  and prompt reply and suggestions you have made

@tomleslie  thanks  for correcting me and  can we collect the real part  and  plot  the graphs