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

Hellow maple users, I am getting an error while solving system of differential equations analytically. Please help to recify the error. Thanks in advance. Here is my codes;

restart:
with(DETools):
# S, N  are constant
Eq1:=diff(u(y),y,y)-u(y)=C(y):
Eq2:=diff(T(y),y,y)=u(y)-diff(u(y),y)^2-u(y)^2+S*T(y)+N*T(y):
Eq3:=diff(C(y),y,y)-C(y)=0:
desys:={Eq1,Eq2,Eq3};ics:={u(0)=0,D(u)(0)=h,T(1)=h,D(T)(0)=0,C(1)=h,D(C)(0)=0}:
combine(dsolve(desys union ics,{u(y),T(y),C(y)}));

Dear maple user how to rectify the error  in solving the coupled differential equation using homotropy perturbation method and direct differentiations and compare the two result by plotting the graphs :

restart:
with(PDEtools):
L:=4:Nb:=1:Nt:=1:#k is some constant
f(x):=sum((p^i)*f[i](x),i=0..L):  
g(x):=sum((p^i)*g[i](x),i=0..L):
HO1:=(1-p)*(diff(f(x),x,x))+p*(((1/x)*(diff(f(x),x))+Nb*((diff(f(x),x))*diff(g(x),x))+Nt(diff(f(x),x)^2))):                                                   
expand(%):                                                                                                                                          
collect(%,p):
HO2:=(1-p)*(diff(g(x),x,x))+p*(((1/x)*(diff(g(x),x))+(Nb/Nt)*((diff(f(x),x,x))+(1/x)*diff(f(x),x)))):                                                                                        
expand(%):                                                                                                               
collect(%,p):                                                                                                              
HO2:=%:                                                                                                                 
declare(f(x),g(x),prime=x): 
for i from 0 to L+1 do equa[1][i]:=coeff(HO1,p,i)=0 end  do:                                     
for i from 0 to L+1 do equa[2][i]:=coeff(HO2,p,i)=0 end  do:           
con[1][0]:=f[0](0)=(h(x)/64),(D(f[0]))(0)=0:                                                                                                                  
con[2][0]:=g[0](0)=-(k-h(x)^2/4),(D(g[0]))(0)=0: 
for j from 1 to L do:                                                                                                    
 con[1][j]:=f[j](0)=h(x),(D(f[j]))(0)=0:                                                                           
 con[2][j]:=g[j](0)=h(x),(D(g[j]))(0)=0:                                                                         
 end do;
for i from 0 to L do;                                                                                                     
dsolve({equa[1][i],con[1][i]},f[i](x));                                                                              
f[i](x):=rhs(%);                                                                                                         
 f[i](x):=evalf(%);                                                                                                          
dsolve({equa[2][i],con[2][i]},g[i](x));                                                                            
 g[i](x):=rhs(%);                                                                                                          
g[i](x):=evalf(%);                                                                                                         
end do; 
for u from 0 to L-1 do:                                                                                                
  f[u](_z1):subs(x=_z1,f[u](x));                                                                                    
 g[u](_z1):subs(x=_z1,g[u](x));                                                                                  
 f[u+1](x):=value(simplify(f[u+1](x)));                                                                         
   f[u+1](x):=simplify(%);                                                                                               
g[u+1](x):=value(simplify(g[u+1](x)));                                                                        
g[u+1](x):=simplify(%);                                                                                                
end do:  
f(x):=evalf(simplify(sum(f[n](x),n=0..L))); 
#### direct calculations                      
#direct solve the two equations in terms of f(x) and g(x) for h(x)=e^x where Nt and Nb #are parameters and its takes some values example 1,1
restart:
with(DETools):
with(plots):
with(IntegrationTools):
Nb:=1:Nt:=1:h(x):=e^x:
Eq1 := (diff(f(x),x,x))+(((1/x)*(diff(f(x),x))+Nb*((diff(f(x),x))*diff(g(x),x))+Nt(diff(f(x),x)^2))):   
Eq2 := (diff(g(x),x,x))+(((1/x)*(diff(g(x),x))+(Nb/Nt)*((diff(f(x),x,x))+(1/x)*diff(f(x),x)))):       

Cd1 := f(0) = h(x), (D(f))(0) = 0:
dsys := {Cd1, Eq1}:
dsol := dsolve(dsys, numeric, output = operator):
#dsol(.1):
plots[odeplot](dsol, [x, diff(f(x), x$1)], 0 .. 5, color = green):
Cd2 := g(0) = h(x), (D(g))(0) = 0:
dsys := {Cd1, Cd2, Eq1, Eq2}:
dsol := dsolve(dsys, numeric, output = operator):
plots[odeplot](dsol, [x, f(x)], 0 .. 5, color = red);
plots[odeplot](dsol, [x,g(x)], 0 .. 5, color = black);
                                                                         

Dear maple user,

I have codes for Differential equations while applying one do and end loop i am able to plot the graph of G(x) while same problem with other way of applying do and end loop i am unable to plot. whats wrong with do and end loop. These are codes available in maple primes . while combining i am unable to plot .

any one resolve it.

Thanks in advance . 

restart:
with(DETools):
with(plots):
with(IntegrationTools):
de0 := {
    (1-p)*(diff(f(x),x,x,x))+p*(diff(f(x),x,x,x)+(1/2)*f(x)*(diff(f(x),x,x))),
    (1-p)*(diff(g(x),x$2))/Pr+p*((diff(g(x),x$2))/Pr+(1/2)*f(x)*(diff(g(x),x)))}:

ibvc0 := {f(0),(D(f))(0),(D(f))(5)-1,g(0)-1,g(5)}:
n:=2:

F := unapply( add(b[k](x)*p^k,k=0..n), x ):
G := unapply( add(c[k](x)*p^k,k=0..n), x ):

de := map( series, eval( de0, {f=F,g=G} ), p=0, n+1 ):

for k from 0 to n do

    if k = 0 then
        ibvc := expand( eval[recurse]( ibvc0, {f=F,g=G,p=0} ) ):
    else
        ibvc := { b[k](0), D(b[k])(0), (D@@2)(b[k])(0), c[k](0), D(c[k])(0) }:
    end if:

    sys := simplify( map( coeff, de, p, k ) ) union ibvc:
    soln := dsolve( sys ):
    
    b[k] := unapply( eval( b[k](x), soln ), x ):
    c[k] := unapply( eval( c[k](x), soln ), x ): 

end do:

'F(x)' = F(x)+O(p^(n+1)):
'G(x)' = G(x)+O(p^(n+1)):

Pr:=1:
plot(eval(G(x), p = 1), x = 0 .. 5, color = blue):
###### Same problem with other  way of do and and end loop unable to plot with G(x)
restart:
with(DETools):
with(plots):
with(IntegrationTools):
Pr:=1:
de1 := (1-p)*(diff(f(x), `$`(x, 3)))+p*(diff(f(x), `$`(x, 3))+(1/2)*f(x)*(diff(f(x), `$`(x, 2))));
de2 := (1-p)*(diff(g(x), `$`(x, 2)))/Pr+p*((diff(g(x), `$`(x, 2)))/Pr+(1/2)*f(x)*(diff(g(x), x)));
ibvc := f(0), (D(f))(0), (D(f))(5)-1, g(0)-1, g(5); n := 2; F := unapply(add(b[k](x)*p^k, k = 0 .. n), x); G := unapply(add(c[k](x)*p^k, k = 0 .. n), x);
DE1 := series(eval(de1, f = F), p = 0, n+1);
DE2 := series(eval(de2, g = G), p = 0, n+1);
CO := map(coeffs, eval([ibvc], f = F), p); CT := map(coeffs, eval([ibvc], g = G), p);

for k from 0 to n do IBVC1 := select(has, C*T, c[k]); slv := dsolve({coeff(DE2, p, k), op(IBVC1)}); c[k] := unapply(rhs(slv), x) end do;
G(x) = G(x)+O(p^(n+1));
plot(eval(G(x), p = 1), x = 0 .. 5);
 

Hellow, 

I am unable to combine the graphs. I have a plot structure as

P1:=plots[odeplot](dsol, [x, F(x)], 0 .. 5, color = red);

P2:=plot(eval(F(x), p = 1), x = 0 .. 5, color = blue);

I want to combine two structure and display in same plot

display (P1,P2);

 thank in advance