ode15s extract solutions for every time step

Hi,

I'm running this system of ode with ode15s, and I would like to visualize one of the parameters ( ptot) at each step of the solution ( t). The timespan of the ode is 0:0.1:100, so I should have 1001 points. However, the .txt files where I'm printing ptot and t only have 44 elements each. Could you please help me find the mistake? I'm trying to plot ptot against t because I'm changing a parameter ( Pin) during the solution and would like to see the effect of this change in time.

Here's the code:

tspan=0:0.1:100;
cond=[51.2112  ; 63.8766  ; 60.7979  ; 49.0010 ;  35.3767 ;  28.5718  ; 33.7930 ;  31.1300 ;  30.6594  ; 29.9741 ;  30.2541  ; 29.6828  ; 28.9798 ]; 

           [t,y, ptot] = ode15s(@fun2_23_08_18_v2,tspan,cond,[])

Ptot = 'ptot_ch.txt';
ttot = 't_ch.txt';
pt_read = importdata(Ptot);
tread=importdata(ttot);
plot(tread,pt_read);

%%% FUNCTION

function [dz,ptot]= fun2_23_08_18_v2(t,y)
pt1=y(1); pt2=y(2); pt3=y(3); pt4=y(4);  pt5=y(5);  pt6=y(6);  pt7=y(7);  pt8=y(8);  pt9=y(9);  pt10=y(10);  pt11=y(11);  pt12=y(12);  pt13=y(13);

    if t<10
      Pin=60;
    else
      Pin=70;
    end
     % -----Constants -----
    N=3.38*10^6; k=2.96*10^-7; alphat=2.6*10^-5; chb=0.2; M=30*10^-9; K=5*10^(-8)*10^-3; H=0.42; S0=0.98; 
  v=3* 7.5*10^-6; 
  r=[0.0119850000000000;0.00958500000000000;0.00764000000000000;0.00604000000000000;0.00473000000000000;0.00366000000000000;0.00400000000000000;0.00575500000000000;0.00726500000000000;0.00889500000000000;0.0107250000000000;0.0128500000000000;0.0153850000000000]; %mm
  L=[1.27076497943190;0.932650928622621;0.544932761536915;0.303082765473283;0.161799106136796;0.155424891414508;0.245072221621871;0.475103125625241;0.273016623935407;0.427646038844292;0.634082325832342;0.846354695529459;0.938696601022114]; %mm
  h=[0.00484000000000000;0.00425000000000000;0.00381000000000000;0.00349000000000000;0.00327000000000000;0.00314000000000000;0.000309000000000000;0.00115000000000000;0.00145000000000000;0.00178000000000000;0.00215000000000000;0.00257000000000000;0.00308000000000000];%mm
  mu=[1.19409872390289e-05;1.12760032450214e-05;1.06583134073916e-05;1.00804835896938e-05;9.56162410894012e-06;9.20633512007761e-06;9.29628357913371e-06;9.96996247072375e-06;1.05291656347798e-05;1.10660983739492e-05;1.16035344790804e-05;1.21594980256614e-05;1.27473361251949e-05]; %mmHg*s
  R= [1.8728    3.1728    4.3411    5.8457    7.8705   20.3059   22.6623   10.9961    2.6277    1.9250    1.4161    0.9612    0.5439]; %[mmHg*s/ml]
  pdrop=[0,6.93,5.87,4.02,2.70,1.82,2.35,2.62,1.27,0.61,0.89,1.31,1.78,2.01];
  n=[1;2;4;8;16;32;64;32;16;8;4;2;1]; 

for i=1:1:14
    if i==1
        pb(i)=Pin;
    else
    pb(i)=pb(i-1)-pdrop(i);
    end
end

diffp=diff(pb)*(-1);
pb_s=[pb(1)+pb(2);pb(2)+pb(3);pb(3)+pb(4);pb(4)+pb(5);pb(5)+pb(6);pb(6)+pb(7);pb(7)+pb(8);pb(8)+pb(9); pb(9)+pb(10);pb(10)+pb(11);pb(11)+pb(12);pb(12)+pb(13);pb(13)+pb(14)];

for z=1:1:14
    S(z)=(((pb(z)^3+150*pb(z))^(-1) *23400)+1)^(-1);
end
deltaS=diff(S)*(-1);

for j=1:1:13

      compl(j)=((3*pi*L(j)*r(j)^3)/(2*Ey*h(j)) )*10^-3; %vessel compliance   (ElBouri2018) [ml/mmHg]   
      Vb(j)=(compl(j)/2)*pb_s(j); %[ml]
      q(j)=diffp(j)/R(j);
      Vt(j)= chb*H*q(j)*deltaS(j)/M;
  end

dpt1=1/(alphat*Vt(1))*( (2*pi*K*r(1)*L(1))/h(1) *(1/2*(pb(1)+pb(2)) -pt1) -M*Vt(1));
dpt2=1/(alphat*Vt(2))*( (2*pi*K*r(2)*L(2))/h(2) *(1/2*(pb(2)+pb(3)) -pt2) -M*Vt(2));
dpt3=1/(alphat*Vt(3))*( (2*pi*K*r(3)*L(3))/h(3) *(1/2*(pb(3)+pb(4)) -pt3) -M*Vt(3));
dpt4=1/(alphat*Vt(4))*( (2*pi*K*r(4)*L(4))/h(4) *(1/2*(pb(4)+pb(5)) -pt4) -M*Vt(4));
dpt5=1/(alphat*Vt(5))*( (2*pi*K*r(5)*L(5))/h(5) *(1/2*(pb(5)+pb(6)) -pt5) -M*Vt(5));
dpt6=1/(alphat*Vt(6))*( (2*pi*K*r(6)*L(6))/h(6) *(1/2*(pb(6)+pb(7)) -pt6) -M*Vt(6));
dpt7=1/(alphat*Vt(7))*( (2*pi*K*r(7)*L(7))/h(7) *(1/2*(pb(7)+pb(8)) -pt7) -M*Vt(7));
dpt8=1/(alphat*Vt(8))*( (2*pi*K*r(8)*L(8))/h(8) *(1/2*(pb(8)+pb(9)) -pt8) -M*Vt(8));
dpt9=1/(alphat*Vt(9))*( (2*pi*K*r(9)*L(9))/h(9) *(1/2*(pb(9)+pb(10)) -pt9) -M*Vt(9));
dpt10=1/(alphat*Vt(10))*( (2*pi*K*r(10)*L(10))/h(10) *(1/2*(pb(10)+pb(11)) -pt10) -M*Vt(10));
dpt11=1/(alphat*Vt(11))*( (2*pi*K*r(11)*L(11))/h(11) *(1/2*(pb(11)+pb(12)) -pt11) -M*Vt(11));
dpt12=1/(alphat*Vt(12))*( (2*pi*K*r(12)*L(12))/h(12) *(1/2*(pb(12)+pb(13)) -pt12) -M*Vt(12));
dpt13=1/(alphat*Vt(13))*( (2*pi*K*r(13)*L(13))/h(13) *(1/2*(pb(13)+pb(14)) -pt13) -M*Vt(13));

pt_tot=[pt1;pt2;pt3;pt4;pt5;pt6;pt7;pt8;pt9;pt10;pt11;pt12;pt13];

for l=1:1:13
    pt_weight(l)=pt_tot(l)*Vt(l);
    Vtt(l)=Vt(l)*n(l);
end
Vt_sum=sum(Vtt);
ptot=sum(1/Vt_sum * (pt_weight));

dz=[dpt1;dpt2;dpt3;dpt4;dpt5;dpt6;dpt7;dpt8;dpt9;dpt10;dpt11;dpt12;dpt13];

   %STORE VALUES
  fileID=fopen('ptot_ch.txt','a');
  fprintf(fileID,' %2.9f\n',ptot);
  fclose(fileID);
  fileID=fopen('t_ch.txt','a');
  fprintf(fileID,' %2.9f\n',t);
  fclose(fileID);
  fileID=fopen('Vt.txt','a');
  fprintf(fileID,' %2.9f\n',Vt);
  fclose(fileID);

Answer 1

[Steven Lord]

There's no mistake. From the description of the tspan input argument on the documentation page for the ode15s function:

"If tspan has more than two elements [t0,t1,t2,...,tf], then the solver returns the solution evaluated at the given points. 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. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step."

To anticipate your next question no, there is no way to force ode15s to step exactly to the list of points in your tspan vector. If it calculates that it can satisfy your tolerances while taking larger (and fewer) steps than you specified it will.

Since you want the value of a variable that's not returned from your function to the ODE solver, the easiest way to do that is probably to modify your function to return that value as the third output argument, call ode15s normally, then call your function using the solution from ode15s as input while specifying the third output. This also avoids needing to worry about rejected steps, where the solver realizes it took too long a step and evaluates your function with a smaller step from the previous point. This appears like the solver stepped backwards, but instead it effectively threw away the longer step and replaced it with the smaller.