Fix step bvo4c

 I send my program m-file 
 attach m.file at email.

I got 21 points for the domain xinit=[0.041 0.042] .when I run the program the domain have been 179 point because I use the function' length(sol.x) ' and result is 179 but I've got 21 point set. I guess that is automated meshing and create 179 mesh .Can you be fixed mesh for the number of pint we see at xinit for example at this program when useing length(sol.x) the result is 21 point ?????

function catalyst_cathod clear all clc P_c=5 ; II=0.05 ; a_ref_4=1e-5 ; C1_ref_4=3.4e-6; D_O2=1.22e-6*(0.4^1.5) ; F= 96487 ; RR=8.314 ; elec=0.57 ; proton=0.07 ; alfa1=2 ; cp1=0.928 ; W1=32 ; cpL=4.2 ; W3 =18 ; antalpy1=-326.36 ; k4_eff= 0.01 ; henry_O2=20.14e4 ;

xinit = [0.041,0.04105,0.0411,0.04115,0.0412,0.04125,0.0413,0.04135,0.0414,0.04145,0.0415,0.04155,0.0416,0.04165,0.0417,0.04175,0.0418,0.04185,0.0419,0.04195,0.042]

solinit =bvpinit(xinit, [II; 200; 1e-5; 1.18; 1.2; 1e-6; 1e-7; 353; 1000]); options=bvpset('AbsTol',1e-8,'RelTol',1e-3);

oft_2=-0.2 for i=1:10 nu=i sol=bvp4c(@anod_H, @bvpbc, solinit,options) ; oft_2=asinh(sol.y(2,end)/sol.y(3,end)*(C1_ref_4/a_ref_4)*0.5)*RR*353/(alfa1*F) end

length(sol.x)

x=linspace(0.041,0.042); u=deval(sol,x);

figure(1) plot(x,u(3,:)); title('concentartion_O2')

figure(2) plot(x,u(1,:)); title('curent_im')

figure(3) plot(x,u(4,:),x,u(5,:)); title('phi s and m') legend('phi-m','phi-s')

figure(4) plot(x,u(6,:)); title('flux-o2')

figure(5) plot(x,u(7,:)); title('flux-water')

figure(6) plot(x,u(8,:)); title('temprature')

figure(7) plot(x,u(2,:)); title('dim/dx')

function dudx = anod_H(x,u) dudx = zeros(9,1);

 dudx(1)= u(2) ;
 dudx(2)=(a_ref_4/C1_ref_4)*((II-u(1))/(4*F*D_O2)*2*sinh(alfa1*F/(u(8)*RR)*(u(5)-u(4)))...
                     +alfa1*F/(u(8)*RR)*u(3)*2*cosh(alfa1*F/(u(8)*RR)*(u(5)-u(4)))*((u(1)-II)/elec+u(1)/proton)) ;
 dudx(3)=(II-u(1))/(4*F*D_O2) ;
 dudx(4)= -u(1)/proton ; 
 dudx(5)= (u(1)-II)/elec ;
 dudx(6)=  u(2)/(4*F) ;
 dudx(7)= -2*u(2)/(4*F) ;
 dudx(8)= u(9) ;
 dudx(9)= ((u(6)*cp1*W1+u(7)*cpL*W3)*u(9)+abs(u(2)/(2*F))*u(8)*antalpy1...
                                     -u(2)*(u(5)-u(4))-u(1)^2/proton)/k4_eff ;
end

 function res = bvpbc(ua,ub)

    PHI=ub(5)-oft_2(1,1);

res = [ ua(1)-II
        ub(1)
        ub(3)-0.065*P_c/henry_O2
        ub(4)-PHI
        ub(5) 
        ub(6)+II/(4*F)
        ub(7)-1e-6
        ua(8)-353.5
        ub(8)-353    ];

 end
end