Sir, There are two coupled equations of f[zeta] and eta[zeta]. From these the equation of f[zeta] is second order differential equation and the equation of eta[zeta] is algebraic equation. I want to plot the graph of eta[zeta] with zeta. So f[zeta] is a function which is to be used in eta[zeta].

Perhaps I have misunderstood some part of your attachement.

If I replace your NDSolve and the plots that follow it with this

f = First[f/.NDSolve[{f''[\[Xi]]-1/f[\[Xi]]^3+1/eo f[\[Xi]] phir[\[Xi]] 10^7==0, f[0]==1, f'[0]==0}, f,{\[Xi],0,1}]];
Plot[f[\[Xi]], {\[Xi], 0, 1}, AxesOrigin->{0,0}, PlotStyle->{Thickness[0.006]}, AxesLabel->{f}, PlotRange->All]
Plot[eta[\[Xi]]-Abs[(\[Xi] Exp[b1[\[Xi]]] a1[\[Xi]] Sqrt[pi/2])/2*Exp[b1[\[Xi]]^2/8] (Erf[Sqrt[2]*\[Xi]+b1[\[Xi]]/
   (2 Sqrt[2])]-Erf[Sqrt[2]*(\[Xi]*(1-\[Beta][\[Xi]])-\[Beta][\[Xi]]*tdash)+b1[\[Xi]]/(2 Sqrt[2])])], {\[Xi], 0, 1},
   AxesOrigin->{0,0}, PlotStyle->{Thickness[0.006]}, AxesLabel->{\[Xi]}, PlotRange->All]

then it seems to find the desired function f and it then uses that f to evaluate your next expression.

I am concerned that in your attachment you seem to want that next expression to equal zero, but from the second plot that does not seem possible.

If I have made some mistake and you can explain what that was then I will see if I can correct it.