ODE 45 step size changes the turning point

I face a problem in ODE 45 that when ever i change my step size which is zspan = [0:0.1:1] in this case, the turning point of the graph changes according to it. If i set my step size at 0.1, the turning point occur at 0.1; step size at 0.5, the turning point occur at 0.5.

It makes me feel like the result that i get is very "rigged" according to the step size indicated.

Is there any way to solve this problem?

function F = Val(z,S)
%A-Fe2O3,B-CH4,C-FeO,D-CO2,E-H2O 
Xsf=S(1);
Xgf=S(2);
Tsf=S(3);
Tgf=S(4);

%Explicit Information

%Mass balance parameter in fuel reactor
R=8.314; %Gas constant-J/(mol.K)
d=0.5; %fuel reactor diameter-m
A=(pi*0.25^2); %Cross-sectional area-m2
pg=0.2411; %gas density-kg/m3
vch4=0.0628; %Vol. flow rate of methane-m3/s
Mwch4=16; %Molecular weight of methane-g/mol
Fgof=1.247; %Initial molar flow rate of fuel(methane)-mol/s
pbulk=2200; %Bulk density of reactive material in reactor-kg/m3
ps=2.81; %Molar bulk density of reactive material-mol/m3
kof=1.31*(10^8); %Pre-exponential factor 
Ea=219900; %Activation energy-J/mol
Cgof=Fgof/vch4; %Initial concentration of fuel(methane)-mol/m3
Tgof=723; %Initial temperature of fuel(methane)-K
n=0.681; %Order of reaction
Fsof=3; %Initial Fe2O3 feed rate-mol/s
kf=kof*exp(-Ea/(R*Tsf));
Cgf=(Cgof*(1-Xgf)*(Tgof))/((1+(2*Xgf))*Tgf);

%Energy balance parameter of fuel reactor
pparticle=4000; %Density of oxygen carrier-kg/m3
dp=75*(10^-6); %mean particle diameter-m
Ap=pi*(dp^2); %Surface area of a particle-m2
Vp=(pi*(dp^3))/6; %Volume of particle-m3
Sm=(pbulk*Ap)/(pparticle*Vp);
Vbar=vch4/A; %Superficial gas velocity-m/s
mug=(0.001596+((3.439*10^-5)*Tgf)+((-8.14*10^-9)*Tgf^2))/1000; %gas viscosity-kg/(m.s)
Nre=(dp*Vbar*pg)/mug;
kg=(-23.35+(0.1698*Tgf)+((1.893*10^-5)*Tgf^2))*(10^-3); %thermal conductivity of gas-W/(m.K)
h=(0.33*(Nre^(1/3))*kg)/dp;
Fsup=0.5; %Molar flow rate of support material-mol/s
CpA=110.9362+(32.04714*(Tsf/1000))+(-9.192333*(Tsf/1000)^2)+(0.901506*(Tsf/1000)^3)+(5.433677/(Tsf/1000)^2); %Heat capacity of Fe2O3-J/(mol.K)
CpB=-0.703029+(108.4773*(Tgf/1000))+(-42.52157*(Tgf/1000)^2)+(5.862788*(Tgf/1000)^3)+(0.678565/(Tgf/1000)^2); %Heat capacity of CH4-J/(mol.K)
CpC=45.7512+(18.78553*(Tsf/1000))+(-5.952201*(Tsf/1000)^2)+(0.852779*(Tsf/1000)^3)+(-0.081265/(Tsf/1000)^2); %Heat capacity of FeO-J/(mol.K)
CpD=24.99735+(55.18696*(Tgf/1000))+(-33.69137*(Tgf/1000)^2)+(7.948387*(Tgf/1000)^3)+(-0.136638/(Tgf/1000)^2); %Heat capacity of CO2-J/(mol.K)
CpE=30.092+(6.832514*(Tgf/1000))+(6.793435*(Tgf/1000)^2)+(-2.53448*(Tgf/1000)^3)+(0.082139/(Tgf/1000)^2); %Heat capacity of H2O-J/(mol.K)
CpF=146.5551+(35.91295*(Tsf/1000))+(-0.183978*(Tsf/1000)^2)+(0.031409*(Tsf/1000)^3)+(-3.659941/(Tsf/1000)^2); %Heat capacity of Support material-J/(mol.K)
HOR1=360642; %Heat of reaction in fuel reactor-J/mol

%Differential Equation
dXgfdz=(ps*A*kf*(Cgf^n))/(4*Fgof);
dXsfdz=(-ps*A*kf*(Cgf^n))/Fsof;
dTgfdz=(A*(Sm*h*(Tsf-Tgf)))/(Fgof*(CpB+Xgf*(2*CpE+CpD+CpB)));
dTsfdz=(A*(Sm*h*(Tsf-Tgf)+((Fgof/A)*(dXgfdz)*HOR1)))/((Fgof*(CpA+(Xsf*(2*CpC-CpA))))+Fsup*CpF);

F=[dXgfdz; dXsfdz; dTgfdz; dTsfdz];

end

To call out the function file:

XsfO = 0.9;
XgfO = 0;
TsfO = 800;
TgfO = 773; 
SO = [XsfO XgfO TsfO TgfO];
zspan = [0:0.1:1];
[z,S] = ode45(@val,zspan,SO);

% Plots

figure(1)
hold on
plot(z,(S(:,1)),'-r','displayname','Xsf')
plot(z,(S(:,2)),'-b','displayname','Xgf')
legend;
ylabel('conversion');
xlabel('Length/m');
hold off

figure(2)
hold on
plot(z,(S(:,3)),'-g','displayname','Tsf')
plot(z,(S(:,4)),'-c','displayname','Tgf')
legend;
ylabel('Temperature(K)');
xlabel('Length/m');
hold off

Answer 1

[Stephan]

This is not a problem of Matlab, but to the misunderstanding how ode solvers in Matlab work. Read here:

https://de.mathworks.com/help/matlab/ref/ode45.html#bu00_4l_sep_shared-tspan

To solve this you should use the following notation:

zspan = [0 1];

instead of

zspan = [0:0.1:1];

since the documentation states:

"...However, the solver does not step precisely to each point specified in tspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points in tspan..."

That means, that if you make the steps to big you get inaccurate turning points, if you make the to small you calculate useless points and waste time on calculations. So you should trust the solver to choose the time steps depending on the problem.