Dynamic response ode15s

Hi,

I would like to get the dynamic response of this system of equations when one of the parameters (Pin) changes in time. Meaning that at baseline conditions Pin=60, but if I increase it at time=20 up to Pin=70, I would like to see how long it takes the system to get to steady state and what the plot looks like. The response I get now is very unusual, could you please help me spot the mistake? Thank you.

tspan=0:0.1:100;
  cond=[%sat
       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 ]; %pt 
  Pin=60; %baseline
           [T,Y] = ode15s(@(t,y)fun1(t,y,Pin),tspan, cond,[]);     

    function dy= fun1(t,y,Pin)
    %baseline Pin=60
    dy=zeros(13,1);

    if t>20
    Pin= 70;
    end

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

    Ey=10^4*0.00750062; 
    v=3* 7.5*10^-6; %kinematic viscosity [mmHg s]
    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];

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

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

%Weighted sum of pt
for l=1:1:13
    pt_weight(l)=pt_tot(l)*Vt(l);
end
Vt_sum=sum(Vt);
ptot=sum(1/Vt_sum * (pt_weight));

dy=[dpt1;dpt2;dpt3;dpt4;dpt5;dpt6;dpt7;dpt8;dpt9;dpt10;dpt11;dpt12;dpt13];
%%STORE VALUES
fileID=fopen('JUMP_ptot_ch.txt','a');
fprintf(fileID,' %2.9f\n',ptot);
fclose(fileID);

This is the code I use for the plot which I will also attach here.

figure(2)
TIME= linspace(0,100,44);
jump = 'JUMP_ptot_ch.txt';
H = importdata(jump);
plot(TIME,H,'Linewidth',2)
title('Dynamic pt change with Pin=40 if TIME>20 else Pin=60')
xlabel('time [s]')
ylabel('weighted partial tissue pressure pt change [mmHg]')