In[]:=
ClearAll["Global`*"]a=1.38*10^10;q=1.6*10^-19;ϵ=8.85*10^-12;λ=4*10^-9;μhat=0.8*10^-20;f=-179;β=4.63*10^-44;α=-80*10^-16;M=70*10^6;eqn1=(M-((α^2)/β))u''[x]-f*P''[x]+(α/β)*μ'[x]==0;eqn2=u''[x]+(a/f)*P[x]-(1/f*q)*μ'[x]==0;eqn3=μ''[x]+(q/ϵ)*P'[x]-(1/λ^2)*μ[x]+(1/λ^2)*μhat==0;sol=NDSolve[{eqn1,eqn2,eqn3,μ'[0]==0,μ[0]==-8*10^-21,P[0]==10,P'[0]==0,u[0]==10^-9,u'[0]==0},{μ,P,u},{x,0,10^-9}]μsol[x_]:=μ[x]/.sol[[1]]Psol[x_]:=P[x]/.sol[[1]]usol[x_]:=u[x]/.sol[[1]]
Out[]=
μInterpolatingFunction,PInterpolatingFunction,uInterpolatingFunction
In[]:=
Plot[{Psol[x]},{x,0,10^-9},PlotLegends->{"P(x)"}]
Out[]=
In[]:=
Plot[{μsol[x]},{x,0,10^-9},PlotLegends->{"μ(x)"}]
Out[]=
In[]:=
Plot[{usol[x]},{x,0,10^-9},PlotLegends->{"u(x)"}]
Out[]=
In[]:=
a1=1.38*10^10;q1=1.6*10^-19;ϵ1=8.85*10^-12;λ1=4*10^-9;μhat1=0.8*10^-20;f1=0;β1=4.63*10^-44;α1=-80*10^-16;M1=70*10^6;eqn4=(M1-((α1^2)/β1))u1''[y]+(α1/β1)*μ1'[y]==0;eqn5=a1*P1[y]-(1/q)*μ1'[y]==0;eqn6=μ1''[y]+(q/ϵ1)*P1'[y]-(1/λ^2)*μ1[y]+(1/λ^2)*μhat1==0;sol=NDSolve[{eqn4,eqn5,eqn6,μ1'[0]==0,μ1[0]==-8*10^-21,P1[0]==10,P1'[0]==0,u1[0]==10^-9,u1'[0]==0},{μ1,P1,u1},{y,0,10^-9}]μsol1[y_]:=μ1[y]/.sol[[1]]Psol1[y_]:=P1[y]/.sol[[1]]usol1[y_]:=u1[y]/.sol[[1]]
Out[]=
μ1InterpolatingFunction,P1InterpolatingFunction,u1InterpolatingFunction
In[]:=
Plot[{Psol[x]},{x,0,10^-9},PlotLegends->{"P(x)"}]
Out[]=