Selecting last value of tspan ode45

I'm trying to extract results only for the last value of the tspan in ode45. However, I end up with 2 results for tfinal... I don't understand why this is happening.

Here's the code:

PS. The ode45 is ran in a for loop because I want to solve the set of diff eqs for 130 different values of ABP (which is a parameter in the eqs).

clear all 
clc

%%=========== ODE45 ==================   
ABP= linspace(40,170,131);

     for i=1:1:length(ABP)  %change in ABP at every loop 
         sol = ode45(@first,[0 100], [0 0 0 0 0.205 97.6 12 49.67],[],ABP(i)); %(function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
     end

%%========= FUNCTION ==================
function dvdt = first(t,v,ABP)

xq= v(1); %variables
xc= v(2);
xm= v(3);
xm1= v(4);
Ca=v(5);
P1= v(6);
V_sa= v(7);
P2 = v(8);

   % -----Constants -----
  R_la= 0.4;
  R_sa_b= 5.03;
  R_sv= 1.32;
  R_lv= 0.56;
  P_v= 6;
  V_la=1;
  V_sa_b= 12;
  P_ic= 10;
  k_ven= 0.186;
  P_v1= -2.25;
  V_vn= 28; 
  tau_q= 20; 
  Pa_co2_b= 40; 
  tau_co2= 40; 
  tau1= 2; 
  tau2= 4;
  tau_g= 1; 
  C_a_p=2.87; 
  C_a_n= 0.164;
  g=0.2; 
  E=0.4; 
  K= 0.15; 
  V0= 0.02; 
  q_b=12.5;
  G_q=3; 
  Pa_co2=40; 
  Ca_b= 0.205;
  eps=1;
  u=0.5; 
  Pa_b=100;
  ka=3.68;

deltaCa_p=2.87;
deltaCa_n=0.164;

% =========== System of diff eq =========================================
dxq= (-xq+G_q*( ( ( (P1-P2)/(R_sv+0.5*(R_sa_b*(V_sa_b/V_sa)^2))) -q_b) /q_b) )/tau_q;
dxc=(-xc +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
dxm=(eps*u-tau2*xm1-xm)/tau1^2;
dxm1=xm1;

sum_xqxcxm=xm+xc-xq; %sum of feedback mechanisms

if t==100 % Last value of tspan!
    disp(sum_xqxcxm) %should get 14 displayed values (because i=14) but I get 28!
end

   if (ABP==40)||(ABP==41) ||(ABP==42) ||(ABP==43) ||(ABP==44) ||(ABP==45) ||(ABP==46) ||(ABP==47)||(ABP==48)||(ABP==49)||(ABP==50)||(ABP==51)||(ABP==52)
      delta=deltaCa_p; %because sum_xqxcxm >0 for ABP=40

elseif (ABP==53)||(ABP==54) %add other values later
     delta=deltaCa_n; %sum<0
 end

dCa=0.5*delta*(1-tanh(2*sum_xqxcxm/delta)^2);
dP1= 1/Ca * ((ABP-P1)/(R_la+0.5*(R_sa_b*(V_sa_b/V_sa)^2)) - (P1-P2)/(R_sv+0.5*(R_sa_b*(V_sa_b/V_sa)^2)) -dCa*(P1-P_ic));
dV_sa= dCa*(P1-P_ic) + Ca*dP1;
dP2=1/(1/(k_ven*(P2-P_ic-P_v1))) *((P1-P2)/(R_sv+0.5*(R_sa_b*(V_sa_b/V_sa)^2)) -(P2-P_v)/R_lv);

dvdt=[dxq;dxc;dxm;dxm1;dCa;dP1;dV_sa;dP2];
%%==== print to file====
Rsa= R_sa_b*(V_sa_b/V_sa)^2;
CBF= (P1-P2)/(Rsa*0.5 + R_sv);
CBF_ch= (CBF-q_b)/q_b *100;

if t==100 %save only final solution
    fileID=fopen('results.txt','a');
     fprintf(fileID,' %-5s %-2s %-5s %-5s %-5s %-5s %-5s %-5s %-5s %-5s\n','xq','xc','xm','xm1','Ca','P1','V_sa','P2', 'CBF','CBFchange'); %how to write this only as a header and not repeat every time?
fprintf(fileID,' %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f\n\n',xq,xc,xm,xm1,Ca,P1,V_sa,P2,CBF,CBF_ch)
fclose(fileID);
end
end

[Jan]

Note: The brute clearing of everything "clear all" is rarely useful, but wastes a lot of time with reloading all functions from the disk.

Where do you "end up with 2 results" wand what do you expect instead? Please explain what "this" is in: "I don't understand why this is happening."

[gorilla3]

if you look in the code, you will see that I put a comment next to disp(sum_xqxcxm). When I display it only 1 value should be displayed, instead I get 2 values.

Answer 1

[Torsten]

Function "first" is called several times for the same t-values.

Better call the function "first" once more after ode45 has finished (i.e. after sol = ode45(...)). Or use OutputFcn.

Best wishes

Torsten.

[gorilla3]

I'm sorry but I really don't understand what you mean by the fact that I'm calling it multiple times. I only defined the function in the beginning.

Could you tell me where the mistake in the code is? I want to display disp(sum_xqxcxm) for every ABP value at each end of the timestep. (i.e. in the for loop ABP has 130 different values for which the ode45 is solved in tspan[0 100]. Hence, I should be getting 130 sum_xqxcxm values for t=100)

[Torsten]

You are not calling "first" multiple times for the same value of t, but ode45 does.

Best wishes

Torsten.

[gorilla3]

So how should the correct code look like for obtaining 130 outputs (i.e. one for every value of ABP) and not more?

I'm contacting you here because I tried numerous ways and I don't know how to fix it

[Torsten]

Make a copy of "first", name it "first2" and change the main program to

ABP= linspace(40,170,131);

       for i=1:1:length(ABP)  %change in ABP at every loop 
           [T,Y] = ode45(@first,[0 100], [0 0 0 0 0.205 97.6 12 49.67],[],ABP(i)); %(function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
           t=T(end);
           y=Y(end,:);
           dy=first2(t,y);
       end

Then make the output in first2.

Or call "first" itself with a negative value for t and only make output in "first" if t<0.

Best wishes

Torsten.

[gorilla3]

Here's a clearer explanation of the problem: https://uk.mathworks.com/matlabcentral/answers/372804-for-loop-ode-45-changing-parameter-values

[gorilla3]

I have tried using your method but I keep getting errors:

ABP= linspace(40,170,131);
t=[0 100];

       for i=1:1:length(ABP)  %change in ABP at every loop        
           [T,Y] = ode45(@first,t, [0 0 0 0 0.205 62.516 12 13.201],[],ABP(i)); %(function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
           t=T(end);
             y=Y(end,:);
             dy=first2(t,y);
       end
  %%========= FUNCTION ==================
  function dvdt = first(t,v,ABP)

xq= v(1); %variables
xc= v(2);
xm= v(3);
xm1= v(4);
Ca=v(5);
P1= v(6);
Vsa= v(7);
P2 = v(8);

     % -----Constants -----
    Rla= 0.4;
    Rsab= 5.03;
    Rsv= 1.32;
    Rlv= 0.56;
    Pv= 6;
    Vla=1;
    Vsab= 12;
    Pic= 10;
    kven= 0.186;
    Pv1= -2.25;
    Vvn= 28; 
    tauq= 20; 
    Pcb= 40; 
    tauc= 40; 
    tau1= 2; 
    tau2= 4;
    taug= 1; 
    g=0.2; 
    E=0.4; 
    K= 0.15; 
    V0= 0.02; 
  %   qb=12.5;
  qb=12.5;
    Gq=3; 
    Pc=40; 
    Ca_b= 0.205;
    eps=1;
    u=0.5; 
    Pab=100;
    ka=3.68;

deltap=2.87;
deltan=0.164;

% =========== System of diff eq =========================================
Rsa= Rsab/((Vsa/Vsab)^2);
q=(P1-P2)/(Rsv+0.5*Rsa);
Cv=1/(kven*(P2-Pic-Pv1));
fn=0.3+3*tanh(Pc/Pcb -1.1);

dxq= 1/tauq *(-xq+Gq*((q-qb)/qb));
dxc=1/tauc*(-xc +fn);
dxm=1/tau2*(-tau1^2*xm1-xm);
dxm1=xm1;

x=xm+xc-xq; %sum of feedback mechanisms
dx=dxm+dxc-dxq;

if x>=0
    delta=deltap;
else
    delta=deltan;
end

dCa=dx* 1/(cosh(2*x/delta))^2;
dP1=1/Ca *((ABP-P1)/(Rla+Rsa*0.5) -q - dCa*(P1-Pic));
dVsa= dCa*(P1-Pic)+Ca*dP1;
dP2=1/Cv *(q-(P2-Pv)/Rlv);

dvdt=[dxq;dxc;dxm;dxm1;dCa;dP1;dVsa;dP2];

end
%%========= FUNCTION 2 ==================
function dy = first2(t,v,ABP)

xq= v(1); %variables
xc= v(2);
xm= v(3);
xm1= v(4);
Ca=v(5);
P1= v(6);
Vsa= v(7);
P2 = v(8);

     % -----Constants -----
    Rla= 0.4;
    Rsab= 5.03;
    Rsv= 1.32;
    Rlv= 0.56;
    Pv= 6;
    Vla=1;
    Vsab= 12;
    Pic= 10;
    kven= 0.186;
    Pv1= -2.25;
    Vvn= 28; 
    tauq= 20; 
    Pcb= 40; 
    tauc= 40; 
    tau1= 2; 
    tau2= 4;
    taug= 1; 
    g=0.2; 
    E=0.4; 
    K= 0.15; 
    V0= 0.02; 
  %   qb=12.5;
  qb=12.5;
    Gq=3; 
    Pc=40; 
    Ca_b= 0.205;
    eps=1;
    u=0.5; 
    Pab=100;
    ka=3.68;

deltap=2.87;
deltan=0.164;

% =========== System of diff eq =========================================
Rsa= Rsab/((Vsa/Vsab)^2);
q=(P1-P2)/(Rsv+0.5*Rsa);
Cv=1/(kven*(P2-Pic-Pv1));
fn=0.3+3*tanh(Pc/Pcb -1.1);

dxq= 1/tauq *(-xq+Gq*((q-qb)/qb));
dxc=1/tauc*(-xc +fn);
dxm=1/tau2*(-tau1^2*xm1-xm);
dxm1=xm1;

x=xm+xc-xq; %sum of feedback mechanisms
dx=dxm+dxc-dxq;

if x>=0
    delta=deltap;
else
    delta=deltan;
end

dCa=dx* 1/(cosh(2*x/delta))^2;
dP1=1/Ca *((ABP-P1)/(Rla+Rsa*0.5) -q - dCa*(P1-Pic));
dVsa= dCa*(P1-Pic)+Ca*dP1;
dP2=1/Cv *(q-(P2-Pv)/Rlv);

dy=[dxq;dxc;dxm;dxm1;dCa;dP1;dVsa;dP2];

qch= (q-qb)/qb *100;

figure(1)
  if t==100
      Vv= (1/kven)*log(P2-Pic-Pv1)+Vvn;
       if t==100 %save only final solution
    fileID=fopen('venousV.txt','a');
fprintf(fileID,' %4.3f\n',Vv);
fclose(fileID);

fileID=fopen('Vsa.txt','a');
fprintf(fileID,' %4.3f\n',Vsa);
fclose(fileID);
       end
  end

I'm making the plots only in first2 now, as suggested.

The error is:
Not enough input arguments.

Error in loop1412>first2 (line 158)
dP1=1/Ca *((ABP-P1)/(Rla+Rsa*0.5) -q - dCa*(P1-Pic));

Error in loop1412 (line 12)
           dy=first2(t,y);

[Torsten]

dy = first2(t,y,ABP(i));

And open your output files in the main program.

And remove the if-statement concerning t=100 in first2 ; they are no longer needed.

Best wishes

Torsten.

[gorilla3]

Wonderful, thanks a lot Torsten!

Answer 2

[Jan]

Do not use the parameter as 5th input, because this is deprecated for over 17 years now. See http://www.mathworks.com/matlabcentral/answers/1971.

A bold guess: You overwrite sol in each iteration. Here is one way to collect them instead:

solC = cell(1, length(ABP));
yIni = [0 0 0 0 0.205 97.6 12 49.67];
for k = 1:length(ABP)
   sol{k} = ode45(@(t, y) first(t, y, ABP(k)), [0 100], yIni);
end

Or maybe you want this:

yR = zeros(length(ABP), 8);
for k = 1:length(ABP)
   [T, Y]   = ode45(@(t, y) first(t, y, ABP(k)), [0 100], yIni);
   yR(k, :) = Y(end, :);
end

[gorilla3]

Thank you, but solC ends up becoming an empty array (each cell contains '[]'). I believe there should be a link between sol and solC?

[gorilla3]

Here's a better explanation of the issue:

https://uk.mathworks.com/matlabcentral/answers/372804-for-loop-ode-45-changing-parameter-values

[Jan]

@gorilla3: As you have found out, "sol" is a typo. Call it "solC" instead. I think, you should be able to fix this by your own, and if the result is provided in a variable called "sol", you should be able to find it there also.