mems.mws

 > restart;

 > ode := diff(w(eta),eta,eta)+(n-1)/eta*diff(w(eta),eta)=1/w(eta)^2;

 > n := 2;

 > L0:=1e-3:

 >

 > pp := dsolve({ode, w(L0)=1, D(w)(L0)=0}, numeric):

 > pp(1);

 > with(plots):

 > odeplot(pp, [eta^2/w(eta)^3, 1/w(eta)], eta=L0..100, numpoints=1300);pic1:=%:

 > odeplot(pp, [eta^2/w(eta)^3, 1/w(eta)], eta=100..4000, numpoints=300, color=blue);pic2:=%:

 >

 > display(pic1,pic2, title="Bif diagram in 2D", labels=["lambda", "u(0)"], plot([4/9,y,y=0..1], color=black, linestyle=dot));

 > n := 1: pp1 := dsolve({ode, w(L0)=1, D(w)(L0)=0}, numeric):

 > odeplot(pp1, [eta^2/w(eta)^3, 1/w(eta)], eta=L0..100, numpoints=1000,  title="Bif diagram in 1D", labels=["lambda", "u(0)"]);

 > pic:=NULL: b:=0.1: to 20 do odeplot(pp, [eta/b, eval(1/w(eta), pp(b))*w(eta)], eta=L0..b, numpoints=300, title=sprintf("lambda=%a", b^2*eval(1/w(eta)^3, pp(b)))); pic:=pic,%: b:=b*2: end:

 > display(pic, insequence=true);

 > plot(eta^(2/3), eta=0..1, color=black, thickness=1);

 > display(%,pic): display(%, title="");

 >

 >