Matlab ode15s change parameter value at specific time during solution

I'm trying to change the value of my variable Pin at specific points in time during the ode15s solution, in order to evaluate the dynamic response. But I get the error:

Error using odearguments (line 83)
The last entry in tspan must be different from the first entry.

I believe the error is somewhere here:

t_start=0;
    t=t_start;
    y=cond;

      while idx_seg< length(t_seg)
          idx_seg=idx_seg+1;
          t_end=t_seg(idx_seg);    
      [t_sol,y_sol]=ode15s(@(t,y)f1_v1,[t_start, t_end],y(end,:));
      t = [t; t_sol(2 : end)];
      y = [y; y_sol(2 : end, :)];
      t_start = t_end;
      end

Here's the full code:

      function [t,y]=f2_v4_21_08_18(cond, t_seg, Pin)

           % -----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; 
      bp=0.8; ro=1040*10^6; 
      Ey=10^4*0.00750062; 
      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]; 
      %%%%%

      idx_seg=0;

      function dy= f1_v1(t,y)
      dy=zeros(13,1);

      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);

       % -----Constants -----
  ...
      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];

        for i=1:1:14
          if i==1
              pb(i)=Pin(idx_seg);
          else
          pb(i)=pb(i-1)-pdrop(i);
          end
          pb(i)=pb;
        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
          kin=v;
          compl(j)=((3*pi*L(j)*r(j)^3)/(2*Ey*h(j)) )*10^-3; 
          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

      %%differential eq
      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];

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

      t_start=0;
      t=t_start;
      y=cond;

      while idx_seg< length(t_seg)
          idx_seg=idx_seg+1;
          t_end=t_seg(idx_seg);    
      [t_sol,y_sol]=ode15s(@(t,y)f1_v1,[t_start, t_end],y(end,:));
      t = [t; t_sol(2 : end)];
      y = [y; y_sol(2 : end, :)];
      t_start = t_end;
      end

      for i=1:1:14
          if i==1
              pb=Pin(idx_seg);
          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
          kin=v;
          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

      for m=1:1:13
                    pt_weight(m)=y_sol(m)*Vt(m);
                 end
                 Vttot=Vt.*n;
      Vt_sum=sum(Vttot);
      ptot=sum(1/Vt_sum * (pt_weight));

      end

These are the initial conditions and times I run it with:

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] = f2_v4_21_08_18(cond, [0, 10, 100], [60, 70, 60]);
plot(t, y);

Thanks in advance!

Answer 1

[Torsten]

t_start = 0 und t_seg(1) = 0.

In the first call to ODE15S, you consequently try to integrate from t_start=0 to t_end=0 which is not accepted by the integrator.

Best wishes

Torsten.

[gorilla3]

Thanks! I now get the following error:

Index exceeds matrix dimensions.

Error in f2_v4_21_08_18/f1_v1 (line 32)
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);

I don't know why since I have 13 initial conditions and then 13 equations... Could you please clarify this?

[Torsten]

cond is 13x1, thus y is 13x1 at the start. Consequently y(end,:) is a scalar, not a vector of length 13.