Dear maple user i am facing difficulty to plot the graph for different values of parameter M=2,4 and fixing t=j=0 to 2 and y=i=0 to 4 on x axis and U on y axis. I am unable to plot 2D . I am enclosing the codes and sample graphs.

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;

T[0, j] := j*lambda1;

T[N+1, j] := 0;

C[0, j] := j*lambda1;

C[N+1, j] := 0

end do:

#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]) ;

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);

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:

#

# Determine the unknowns in the system

#

`union`( seq(seq( indets( eq1[i,j], name), i=1..N), j=0..N),

seq(seq( indets( eq2[i,j], name), i=1..N), j=0..N),

seq(seq( indets( eq3[i,j], name), i=1..N), j=0..N)

);

#

# And how many unknowns

#

numelems(%);

#

# And the number of equations

#

numelems(eq1)+numelems(eq2)+numelems(eq3);

#

# So if one supplies values for 'Gr' and 'M', and

# (assuming the equations are consistent), one ought

# to be able to get values for C[1..9, 1..10],

# T[1..9,1..10], and U[1..9,1..10]

#

# As an example below, choos Gr=1.0 and M=2, then the

# following obtains a 'solution` afer a minute or so

#

fsolve( eval( [ seq(seq(eq1[i,j], i=1..N),j=0..N),

seq(seq(eq2[i,j], i=1..N),j=0..N),

seq(seq(eq3[i,j], i=1..N),j=0..N)

],

[Gr=1.0, M=2]

)

);