Question: How to generate the graphs for the following codes

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