@Carl Love 

Sir please give any suggestions how to proceed or any other method can We apply

@Carl Love 

Then sir how to handle if 3 variables X,R , t are there,

@Rouben Rostamian  

thank for pointing the error .Here i have  make some correction as per your comment but still error. I am attaching the file.

system_of_pdes.mw

@Carl Love

Sir per your comment i am trying to solve the system of difference scheme which is given in an article and try to plot the graphs , but i am unable to do so . for your reference i am attaching an article i followed the same schme which is given in article.

 help_article.pdf

@Rouben Rostamian  Thanks

what is the value of L/g

@mmcdara 

Sir from the output of the codes please help to collect U,T and C values.

@tomleslie , @ Preben Alsholm 12698 

Good day have a look at this codes still i am getting error.  

still_error.mw

@tomleslie 

Dear please see this equations

@tomleslie 

Dear sir  please see the  original equations:  Two coupled equations equations are there with bcs. where Nt, Nb are parameters it will take some constant values.these equaions i want to solve .

Eq1 := diff(theta(r,z),r,r))+(((1/r)*(diff(theta(r,z),r))+Nb*((diff(theta(r,z),r))*diff(sigma(r,z),r))+Nt*(diff(theta(r,z),r)^2)));
Eq2 := (diff(sigma(r,z),r,r))+(((1/r)*(diff(sigma(r,z),r))+(Nb/Nt)*((diff(theta(r,z),r,r))+(1/r)*diff(theta(r,z),r))));
Cd1 := theta(r,z)(h(z),0) = 0, (D(theta(r,z))) = 0;
Cd2 := sigma(r,z)(h(z),0) = 0, (D(sigma(r,z))) = 0;

@tomleslie  Thanks sir for nice explanation. In this way we have to write the boundary conditions sir

con[1][0]:=f[0](h)=0,(D(f[0]))(0)=0:                                                                                                                  
con[2][0]:=g[0](h)=0,(D(g[0]))(0)=0: 
for j from 1 to L do:                                                                                                    
con[1][j]:=f[j](h)=0,(D(f[j]))(0)=0:                                                                           
con[2][j]:=g[j](h)=0,(D(g[j]))(0)=0:

Directely i am getting an error. Please see the  codes

with(DETools):
with(plots):
with(IntegrationTools):
Nb:=1:Nt:=1:h(x):=exp(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)))):       

# Boundary condition f(x)=0, g(x)=0 at x=h(x)" where h(x) is some function of x", 

 #f '(x)=0, g'(x)=0 at x=0

Cd1 := f(h(x)) = 0, (D(f))(0) = 0 :
dsys := {Cd1, Eq1}:
dsol := dsolve(dsys, numeric, output = operator):
plots[odeplot](dsol, [x, diff(f(x), x$1)], 0 .. 5, color = green):
Cd2 := g(h(x)) = 0, (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);

Thanks 

@tomleslie 

Thanks for clear and patience  explanation.

My main intention is to solve coupled differentiation equations solve directly and the same equations i want to solve by using Homotropy perturbation method and compare the two results. By one method of HPM i can produce the graph with more error.

By another method of  HPM codes i can produce the graph of f(x) with less error but g(x) i can't produce the graphs. As i have seen the article by D.D Ganji " Application of He's homotropy perturbation method to boundary layer flow and convection heat transfer over a flat plate" 2007. where he has solve  and compare the results.

 Hope you  can solve my query.

Test3.mw

@tomleslie 

Sir same equations with differerent do else why we can not calculate G(x) . In one loop we can figure out and in other we cant.  I want to write this codes 

.........................................................

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

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

ibvc := f(0), (D(f))(0), (D(f))(5)-1, g(0)-1, g(5)

for k from 0 to 1 do IBVC2 := select(has, CT, c[k]); slv1 := dsolve({coeff(DE2, p, k), op(IBVC2)}); c[k] := unapply(rhs(slv1), x) end do;
G(x) = G(x)+O(p^(n+1));
   I above codes working fine and produces the graph.

Test_2.mw

@tomleslie 

Dear sir,

I want to apply like this conditions

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

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:

@tomleslie  Dear sir please see the files I have same equations  with one method i am getting a solution where as with different method i am unable to plot, i.e I can able to plot for F (x), G(x)  with one method and not able to  produce G(x) with other.

I am uploading the file

thanks in advance

Test_1.mw