5.79 How to numerically solve a BVP ode and plot the solution?

This one has one solution

eq:=diff(u(z),z$2)+(k-1)*diff(u(z),z)/z+lambda*exp(u(z))=0; 
sol:=dsolve({subs({k=1,lambda=2},eq),u(0)=1,u(1)=0},numeric,u(z), 
            method=bvp[midrich],'abserr'=0.001); 
plots[odeplot](sol);
 

This solved coupled ODE’s, so there are 2 solutions. Say \(x_1(t)\) and \(x_2(r)\), It is a little tricky to plot all solutions generated, but here is an example

restart; 
R := 0.4; px := 32000; Mm := 0.1; Ds := 9; DO2 := 7.2; YXS := 0.3; KS := 10; 
Sp := 30; Cb := 8; KO2 := 0.2; R0 := 0.000001; YXO := 0.42857; 
Vs := px*1/YXS*(Mm*x2(r))/(KS + x2(r))*x1(r)/(KO2 + x1(r)); 
Vo := px*1/YXO*(Mm*x2(r))/(KS + x2(r))*x1(r)/(KO2 + x1(r)); 
 
eqs := diff(x1(r),r$2) + 2/r*diff(x1(r),r)= Vo/DO2, 
diff(x2(r),r$2) + 2/r* diff(x2(r),r)= Vs/Ds; 
ic:=D(x1)(R0)=0,x1(R) = Cb,D(x2)(R0)= 0, x2(R) = Sp; 
sol:=dsolve({eqs,ic},numeric,{x1(r),x2(r)},'abserr'=.52,'maxmesh'=1000,output=listprocedure);
 

And now to plot do

x1Sol:=rhs(sol[2]); 
plot(x1Sol(r),r=0..0.4); 
 
x2Sol:=rhs(sol[4]); 
plot(x2Sol(r),r=0..0.4);