@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(`δ1`), typeset(conjugate(`Δ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

@mmcdara 

these are the complete equations given in the article

eq:= '1/r*diff(r*diff((1-(mu/eta)*alpha^2*1/r*diff(r*diff(u(r), r), r)), r), r)'='diff(P(z), z) ':
#boundary conditions:
'diff(u(r), r)=0':  # at r=0
'diff(u(r), r,r)-(mu/eta)*diff(u(r), r)=0': #AT r=h(z)
'diff(u(r), r,r)-(mu/eta)*diff(u(r), r)': # is finite AT r=0
'u':# is finite at r=0
'u':# is finite at r=h(z)
 

@Mariusz Iwaniuk 

thanks for providing the article

@mmcdara  in the article they got the results in terms of bessels

@mmcdara   thanks 

yes u is a function of r and z, whre P is dp/dz and h(z) =1-3*delta*(11*(z-d)), where delta and d is constants , r  is a function of z.

this equations come from the article " https://www.banglajol.info/index.php/JNAME/article/view/31165"

i want to solve and plot it completely using maple,

i