# Question:How to generate the graphs for the following codes

## Question:How to generate the graphs for the following codes

Maple

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

﻿