@Carl Love  thanks for  responding the comment i further trying to understand the equations  Eqs 23-31. further i dug  one more article of gupta el. al there they have given some  methodology . As far as my knowledge i think  they have used matlab to solve those equations .  As   ( Bagh Ali 5 ) particulary want procedure to solve FEM using maple by  variational formula . so it is very helpful that  a procedure can be written on maple so that comparision with analytical , FEM  and error also can be found with the help of graph and numerically.

FEM_paper2.pdfFEM_paper2.pdf

@Carl Love  sorry for reminding again and again , i am also enthusiastic about how to solve non- linear differential equations by  using FEM( maple codes). what ever codes you have applied to solve the non-linear differential equation i have applied the same codes to some other problem still i am getting error . I am attaching the complete paper which i have try to solve , but particularly  ( Bagh Ali 5 ) in his comment ask to solve the non-linear differential  equation  using FEM by maple.

I am also want to interested to see how the nonlinear ODE can be solved using  FEM in maple.

here the maple codes and paper attaching

restart:
Digits:= 15:

ODEs:= [
   #Eq 1/2/3:
   (1+K)*diff(u(y),y$2) + K*diff(N(y),y) + theta(y)-M=0,
   #Eq 10/19/25:
   (1+K/2)*diff(N(y),y$2) -K*(2*N(y)+diff(u(y),y$1))  =  0,
   #Eq 11/20/26:
   diff(theta(y),y$2)=  0
]:                                                            BCs:= [
   #Eq 4:
   u(0) = 0,theta(0) = R, N(0) =0, 
   #Eq 5:
   u(1) = 0,theta(1) = R, N(1) = 0 
   
];                                                               param_names:= [K,M, R];                                       <ODEs[]>;

     BCs:=[uApplyFunction(0)=0,?ApplyFunction(0)=R,NApplyFunction(0)=0,u

       ApplyFunction(1)=0,?ApplyFunction(1)=R,NApplyFunction(1)=0]
                            param_names:=[K,M,R]
                           Matrix(%id = 163720632)
Solve:= 
   subs(
      _P= param_names,
      proc({
         K::realcons:= 5,    
         R::realcons:= 0.5,    
         M::realcons:= 5
      })
         userinfo(1, Solve, param_names=~ _P);
         dsolve(
            eval(
               convert(ODEs, set) union convert(BCs, set),
               param_names=~ _P
            ),
            numeric
         )
      end proc
   )
;
proc({K::realcons := 5, M::realcons := 5, R::realcons := 0.5})  ...  end;
infolevel[Solve]:= 1:

P:= M:
vals:= [0, 5, 10]:
sols:= [seq(Solve(P= v), v= vals)]:
colors:= [red, green, blue, black,orange]:
for F in [u,N](y) do
   print(plots:-display(
      [seq(
         plots:-odeplot(sols[k], [y,F], color= colors[k], legend= [P= vals[k]]),
         k= 1..nops(vals)
      )],
      labeldirections= [horizontal,vertical]
   ))
od:
Solve: [K = 5 M = 0 R = .5]
Solve: [K = 5 M = 5 R = .5]
Solve: [K = 5 M = 10 R = .5]

sols(1):
P:= K:
vals:= [0, 5, 10]:
sols:= [seq(Solve(P= v), v= vals)]:
colors:= [red, green, blue, black,orange]:
for F in [u,N](y) do
   print(plots:-display(
      [seq(
         plots:-odeplot(sols[k], [y,F], color= colors[k], legend= [P= vals[k]]),
         k= 1..nops(vals)
      )],
      labeldirections= [horizontal,vertical]
   ))
od:
paper.pdf

@waseem 

Any idea how to apply FEM to solve

@Carl Love

Did you got any clue  to write the codes in maple

@Carl Love 

Any idea ,   can you give some reference in maple to write the codes for FEM  for second ODE 

@Carl Love complete_paper.pdf

please see the whole paper, it is better to use the same technique for different elements , if maple can solve using the technique mentioned in the paper, so number of bvp can be solved.

@Carl Love pages.pdf

@mmcdara  thanks , can we plot contour plot by fixing Nb=0.1, where delta2 varies from 0.02..0.1 . 

for different value of Gamma2=10,20,30,40,50,-10.

For Nb=0.1, Gamma2=10 one contour plot, similarly

                     Gamma2=20  another contour plot    so on. total 6 contour plot

here is the codes

restart;
  h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
  K1:=((4/h(z)^4)-(sin(alpha)/Gamma2)-h(z)^2+Nb*h(z)^4):
  lambda:=unapply(Int(K1,z=0..1), Gamma2):
  L1:=0.2:
  d1:=0.2:
  alpha:=Pi/6:
  with(plots):
  display
  ( Matrix(
    (3,2, [ seq
        ( plot
          ( [ seq
              ( eval(lambda(Gamma2), Nb=j),
                j in [0.1,0.2,0.3]
              )
            ],
            delta2=0.02..0.1,
            legend=[Nb=0.1,Nb=0.2,Nb=0.3],
            labels=[typeset(`&delta;1`), typeset(conjugate(`&Delta;p`))],
            title=typeset("Effect of ", ''alpha'', " when ", Gamma,"2=", Gamma2)
          ),
          Gamma2 in [10,20,30,40,50,-10]
        )
      ]
    ))
  );
 

@tomleslie 

excellent, thanks a lot

@tomleslie 

sir if we edit the last line we will get the graphs like this i am attaching the sample graph. I want the title  in the graph sample like that and on y axis i want delta p (bar)

*** Please see the title and vertical line label

plots:-display(Vector[row]([seq(plot([seq(eval(lambda(Gamma2), Nb=j), j=[0.1,0.2,0.3])], delta2=0.02..0.1, legend=[Nb=0.1,Nb=0.2,Nb=0.3],labels = ["d1 ", "D"],labelfont = ["SYMBOL", 12], title=( 'Gamma2'=Gamma2)), Gamma2=[10,20,30])]),axes=boxed);

@Kitonum thanks 

@Rouben Rostamian  

can we get the solution of Differential equation with out using boundary condition  in terms of 

Theta(r,z)=,,,,, and Sigma(r,z)

@Rouben Rostamian  

Some body solve this by applying homotropy perturbation method and they got the result as 

Theta(r,z)=(r^2-h^2)/4

Sigma(r,z)=-(r^2-h^2)/4

Note: can we solve this by series also, please give suggestions

@Rouben Rostamian  

 Suppose if we did not  apply boundary condition then what is the solution

@mmcdara Thanks a lot for your complete explanation.

if we treat dp/dz = P=constant  than can we calculate the value of u.

i will be thankful to you if you help  me to  write the  code how to calculate the value of u by fixing P=alpha=eta=d=mu= some constant  and h(z) =1-3*delta*(11*(z-d)), i want to plot a graph by taking u on vertical axis, delta on horizonal axis by taking delta=0...0.01, where z=0 to 1,by changing mu=0,0.01,0.02. sample graph