How to call multiple functions in ode 45?

I'm trying to solve 2 systems of differential equations using ode45, however I get the error code:

Error using odearguments (line 21) When the first argument to ode45 is a function handle, the tspan argument must have at least two elements.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in CBF_2006_f_v1 (line 49)
    [ABP,P1,P2,xq,xm,xc,xm1,Cav] = ode45(@vsa, @aut, tspan, var0);

There are 2 systems: 4 equations each. The unknowns are (V_sa,P1,P2,xq,xm,xc,xm1,Cav), all the rest are parameters that have been declared at the beginning of the code.

The code main code that recalls the functions is:

%Parameter declaration .... (Rsa,V_sa_b,P_ic ...)

%%Diff eq system

      ABP= linspace(40,170,131);

      delta=deltaCa_p;

      for i=1:1:length(ABP)
      tspan=[0 100];

      var0=[12; 97.6; 49.67; 0; 0; 0; 0; 0.205];
      [ABP,P1,P2,xq,xm,xc,xm1,Cav] = ode45(@vsa, @aut, tspan, var0);

      end

The 2 functions are stored in separate files and their code is:

function dvdt = vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1)

dvdt = [(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1)))*(var(2)-P_ic) + var(8)*(( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8)) ;
    ( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8);
                  ((var(2)-var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) +R_sv)-(var(3)-P_v)/R_lv)*1/(1/(k_ven*(var(3)-P_ic-P_v1))) ;
                      0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))] ;

V_sa= var(:,1);
P1= var(:,2);
P2 = var(:,3);
Cav=var(:,8);

xq= var(:,4);
xc= var(:,5);
xm= var(:,6);
xm1= var(:,7);

dvdt=[V_sa;P1;P2;Cav; xq;xc;xm;xm1;Cav];
end

and

    function dxdt= aut(ABP,var,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)

    dxdt=[                                     (-var(4)+G_q*( ( (var(2)- var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) + R_sv) ) -q_b /q_b) )/tau_q ;
                                               (-var(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
                                                                  var(7);
                                                                     (eps*u-tau2*var(7)-var(6))/tau1^2 ];
    V_sa= var(1);
    P1= var(2);
    P2 = var(3);
    xq= var(4);
    xc= var(5);
    xm= var(6);
    xm1= var(7);

    dxdt=[P1;P2;xq;xc;xm;xm1];
    end

Could you please help me find the mistakes?

Answer 1

[Star Strider]

I would nest them into one function:

bothfcns = @(ABP, var) [@(ABP,var) vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1); @(ABP,var) aut(ABP,var,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)];

Then use ‘bothfcns’ as the function argument to ode45.

NOTE — I cannot run your code to test this, so I am listing this as UNTESTED CODE.

[gorilla3]

Thank you very much! No error appears now, however it does not show me the solutions of the equations anywhere. (i.e. the values of the variables V_sa,P1,P2,Cav,xq,xm,xc,xm1). It's not clear to me why... Could you help me further?

[Star Strider]

My pleasure.

The ODE solvers return only the independent variable (as a column vector), and the matrix of integrated dependent variables as a function of the independent variable, each in a column of that matrix. If you want specific values from them, you either have to identify them from the matrix, or calculate them from those values.

That is likely the only other help I can provide, since I do not know what you are doing in your differential equations.

[gorilla3]

Thanks again for this elaborate response.

Alright, I see there's a "bothfncs" element in the workspace, however when I try opening it it says only 1x1 function_handle. What I need, is extracting the results of xq,xc,xm at every time step and summing them. (ie. x_res=xm_res+xc_res+xq_res) Also, at every time step the result for each should be a vector so I would only need to consider the last element of the vector.

[Star Strider]

My pleasure.

I created ‘bothfcns’ to vertically concatenate ‘vsa’ and ‘aut’ so they would work in ode45. It is only a function handle, and it is not possible to get anything else from it, especially the result of the ode45 integration. Those results are in the independent variable ‘ABP’ and the dependent variable matrix ‘var’.

The way I read your code:

xq=var(:,12);
xc=var(:,13);
xm=var(:,14);

(you need to verify that the column references are correct) are in the integrated output after ode45 completes.

I have absolutely no idea what ‘xm_res’ and the rest are, since they appear in your last comment for the first time. However to sum the corresponding columns in ‘var’, just select the columns and sum across the rows:

sum_xqxcxm = sum(var(:, 12:14), 2);

(again assuming I have the columns correct). The result, ‘sum_xqxcxm’ will be a column vector. The last element is simply ‘sum_xqxcxm(end)’. Alternatively, if you only want to sum the last element of each vector, that would be:

sum_xqxcxm = sum(var(end, 12:14), 2);

It is neither appropriate nor likely possible to do the column summation while ode45 is calculating them. You have to wait until it is finished, and then do the calculations on the results.

[gorilla3]

Hi,

this makes a lot of sense, thanks! However, I don't understand how you got to the numbers 12-13-14, because counting them in "bothfncs" I get 6-7-8.

I tried running the code with sum_xqxcxm but I get an Error: "not enough input arguments". Could you please help?

[Star Strider]

My pleasure.

Go with 6:8. You know your code better than I do.

This works for me:

var = randi(9, 5, 18);                          % Create ‘var’
sum_xqxcxm = sum(var(:, 6:8), 2);

You did not post the code you used to calculate ‘sum_xqxcxm’, so I cannot tell you why it failed.

[gorilla3]

Hi,

this is the code I used:

sum_xqxcxm = sum(var(end, 12:14), 2);

Why would var be a set of random numbers?it should be the vector of solutions of the ode45...

Also, as you suggested, I inserted 'bothfncs' as the function argument to ode45 and this is the error message:

Error using feval
Undefined function or variable 'bothfncs'.

For completion here's the full code:

%%Parameters
R_la= 0.4; % mmHg s/ml
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; %[/ml] elastance coeff for venous compliance
P_v1= -2.25; % pressure offset for venous compl
V_vn= 28; %offset of vG_q= 3; %%gain of flow based feedback mechanism
tau_q= 20; %[s] time const of flow based feedback mechanism
Pa_co2_b= 40; % baseline arterial CO2 pressure
tau_co2= 40; %time const for arterial co2 p
tau1= 2; %1st time const of neural feedback
tau2= 4;
tau_g= 1; % t const for tissue O2 concentration
C_a_p=2.87; %ampl of + change in art comp [ml/mmHg]
C_a_n= 0.164;
g=0.2; %baseline non dim tissue O2 concentration
E=0.4; %baseline O2 extraction fraction
K= 0.15; %scaling between activation and O2 demand
V0= 0.02; %scaling factor for BOLD response
q_b=12.5;
G_q=3; %gain of flow based feedback mechanism [ml/mmHg]
Pa_co2=40; 
Ca_b= 0.205;
eps=1;
u=0.5; %0 for steady state
Pa_b=100;
ka=3.68;
deltaCa_p=2.87;
deltaCa_n=0.164;
%%Diff eq system
     ABP= linspace(40,170,131);  

      delta=deltaCa_p;

      tspan=[0 100];

      var0=[12; 97.6; 49.67; 0; 0; 0; 0; 0.205];

      bothfcns = @(ABP, var) [@(ABP,var) vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1); @(ABP,var) aut(ABP,var,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)];
  [ABP,var] = ode45('bothfncs',tspan, var0);

sum_xqxcxm = sum(var(:, 12:14), 2);

function dvdt = vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1)

      for i=1:1:length(ABP)
  dvdt = [(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1)))*(var(2)-P_ic) + var(8)*(( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8)) ;
      ( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8);
                    ((var(2)-var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) +R_sv)-(var(3)-P_v)/R_lv)*1/(1/(k_ven*(var(3)-P_ic-P_v1))) ;
                        0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))] ;

V_sa= var(:,1);
P1= var(:,2);
P2 = var(:,3);
Cav=var(:,8);

xq= var(:,4);
xc= var(:,5);
xm= var(:,6);
xm1= var(:,7);

dvdt=[V_sa;P1;P2;Cav; xq;xc;xm;xm1;Cav];
end
end

function dxdt= aut(ABP,V_sa,P1,P2,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)

dxdt=[                                     (-var(4)+G_q*( ( (var(2)- var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) + R_sv) ) -q_b /q_b) )/tau_q ;
                                           (-var(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
                                                              var(7);
                                                                 (eps*u-tau2*var(7)-var(6))/tau1^2 ];
V_sa= var(1);
P1= var(2);
P2 = var(3);
xq= var(4);
xc= var(5);
xm= var(6);
xm1= var(7);

dxdt=[P1;P2;xq;xc;xm;xm1];
end

[Star Strider]

‘Why would var be a set of random numbers?’

In your code, it wouldn’t. I created it to test my code to be certain I hadn’t made a typographical error.

Undefined function or variable 'bothfncs'.

Please check your spelling:

bothfcns = @(ABP, var) [ ...

[ABP,var] = ode45('bothfncs',tspan, var0);

should be:

[ABP,var] = ode45(bothfcns,tspan, var0);

[gorilla3]

I already tried removing the '' and I still get an error:

Undefined function or variable 'bothfncs'.

Error in CBF_2006_f_v4 (line 48)
[ABP,var] = ode45(bothfncs,tspan, var0);

[Star Strider]

You have to spell the function name correctly in the ode45 call.

Try this:

[ABP,var] = ode45(bothfcns,tspan, var0);

[gorilla3]

Huh, I'm very sorry for that...still:

Error using vertcat
Nonscalar arrays of function handles are not allowed; use cell arrays instead.

Could you please be so kind and look over the full code I sent you above? It would really mean a lot. I was asked to specifically use this "function method" in ode45 and I'm struggling since it's my first time.

[Star Strider]

I managed to get ‘vsa’ working by using single subscripts only to refer to elements of ‘var’.

Note that both ‘dvdt’ and ‘dxdt’ return (4x1) column vectors.

In ‘aut’, you use several variables you apparently have not defined, so even when I add them to the argument list, the function fails with a ‘Not enough input arguments’ error. You also did not include ‘var’ as the second argument to it, as required, Adding that throws other errors when I tried to run it independently.

As an initial approach, get each of your two differential equations working with ode45 first. When they run, try them with my original ‘bothfcns’ anonymous function.

Your code is so complicated that I can’t even begin to figure it out. I will help you with it as I can. You have to find out what problems you are having with it.

Also, my original function vertical concatenation idea works correctly in this test code:

f1 = @(x) [x(1).^2; x(2).^3];
f2 = @(x) [sqrt(x(1)); x(1)*x(2)];
f3 = @(x) [f1(x); f2(x)];

x = [2; 3];
y = f3(x);

I have no idea what the problem was when you were unable to use my ‘bothfcns’ anonymous function. I cannot test it with them because I cannot get ‘aut’ to work.

[gorilla3]

Thank you for all your advice :) I'll try harder to solve it all.

[Star Strider]

As always, my pleasure.

I will help as I can. It would help if I knew what you are doing and what you want to simulate.

[Greg Heath]

Why are you using the name of the MATLAB function var?

 help var
 doc var

?

Hope this helps.

Greg

[gorilla3]

I would like to solve the 2 systems of equations, then (from aut) extrapolate the solutions for xc,xm,xq and sum them. Subsequently I would like to check whether the sum is positive or negative and according to the sign, allocate a value to the parameter delta. If the sum is + then delta=2.8, else delta=0.16. Delta is used in system vsa.

I fixed the function aut, and re-named it "second"(because I read file name and function name should not be the same), so it's working now:

function dxdt= second(ABP,var, V_sa,P1,P2,Cav,xq,xm,xc,xm1,R_sa_b,V_sa_b,R_sv,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)

dxdt=[                                     (-var(4)+G_q*( ( (var(2)- var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) + R_sv) ) -q_b /q_b) )/tau_q ;
                                           (-var(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
                                                              var(7);
                                                                 (eps*u-tau2*var(7)-var(6))/tau1^2 ];
V_sa= var(1);
P1= var(2);
P2 = var(3);
xq= var(4);
xc= var(5);
xm= var(6);
xm1= var(7);

dxdt=[V_sa;P1;P2;xq;xc;xm;xm1];
end

__________________

Thank you Greg, but that's not the issue.

[Star Strider]

The file name for a function should be the same as the function name.

I’m still not certain that I understand what you’re doing.

As I read it, ‘aut’ calculates as ‘dxdt’ a (4x1) column vector, so it should have a vector of 4 initial conditions and return the (4x1) vector initially calculated by ‘dxdt’ only. (You can calculate the other variables from ‘dxdt’ outside the function.) As you wrote it, it overwrites ‘dxdt’ and returns the inputs (the ‘var’ vector) without actually calculating anything with them. This is not going to give you the result you want (that being the integrated differential equation).

[gorilla3]

Hi, I improved and simplified the code. Here it is:

clear all 
clc

%%Diff eq system

function second

[ABP,par] = ode45('first',[0 100], [12 97.6 49.67 0.205 0 0 0 0]);

sum_xqxcxm = sum(par(:, 4:6), 2);

function dvdt = first(ABP,par)
    %%Parameters
R_la= 0.4; % mmHg s/ml
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;
 ABP= linspace(40,170,131); 

       for i=1:1:length(ABP)

% %calculate XSUM first and find 'i' values
% if i== (1||2||n...)
% delta=deltaCa_p; 
% elseif i==(4||5...100)
% delta=deltaCa_n;
% end
delta=2.15 % this is just to make it run but needs to be fixed, see explanation below
dV_sa= dCa*(par(2)-P_ic) + par(8)*dP1;
dP1= 1/par(8) * ((ABP(i)-par(2))/(R_la+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) - (par(2)-par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) -dCa*(par(2)-P_ic));
dP2=1/(1/(k_ven*(par(3)-P_ic-P_v1)))*((par(2)-par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) -(par(3)-P_v)/R_lv);
dCa=0.5*delta*(1- ((cosh(4*(par(6)+par(5)-par(4))/delta))-1)/(cosh(4*(par(6)+par(5)-par(4))) +1));

dxq= (-par(4)+G_q*( ( (par(2)- par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) ) -q_b /q_b) )/tau_q ;
dxc=(-par(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
dxm1=par(7);
dxm=(eps*u-tau2*par(7)-par(6))/tau1^2;

V_sa= par(1);
P1= par(2);
P2 = par(3);
Cav=par(8);

xq= par(4);
xc= par(5);
xm= par(6);
xm1= par(7);
par=[V_sa;P1;P2;Cav;xq;xc;xm;xm1];

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

end

As you can see, at each time step the value of ABP changes (for i=1:1:length(ABP)). I would now like to extract the solutions for xq,xc,xm at each time step t. Then I would like to sum these 3 solutions for each time step and print them in an array.

sum_xqxcxm = sum(par(:, 4:6), 2);

Then check for each value int he array if the value is >0 or <0. Assign delta=2.87 if positive and delta=0.16 if negative. Where delta is a parameter in the system of equations.

[gorilla3]

I figured I could write the part of the delta code like this:

if sum_xqxcxm(i)<0
   delta=deltaCa_n;
elseif sum_xqxcxm(i)>0
   delta=deltaCa_p;      
end

[gorilla3]

Also, when I remove the outermost function second (i.e. just delete its declaration in the beginning and the "end" at the bottom) the code isn't working. Could you tell me why?

[Torsten]

Why do you use the loop

for i=1:1:length(ABP)
...
end

?

The array "dvdt" is only filled with the values for the last loop index i=length(ABP).

Best wishes

Torsten.

[gorilla3]

Because ABP is

ABP= linspace(40,170,131);

Hence, I would like to calculate the differential equations for all the different ABP values.

[gorilla3]

So I guess I should change it to:

dV_sa(i)= dCa*(par(2)-P_ic) + par(8)*dP1;
dP1(i)= 1/par(8) * ((ABP(i)-par(2))/(R_la+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) - (par(2)-par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) -dCa*(par(2)-P_ic));
dP2(i)=1/(1/(k_ven*(par(3)-P_ic-P_v1)))*((par(2)-par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) -(par(3)-P_v)/R_lv);
dCa(i)=0.5*delta*(1- ((cosh(4*(par(6)+par(5)-par(4))/delta))-1)/(cosh(4*(par(6)+par(5)-par(4))) +1));

dxq(i)= (-par(4)+G_q*( ( (par(2)- par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2))) -q_b /q_b) )/tau_q ;
dxc(i)=(-par(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
dxm1(i)=par(7);
dxm(i)=(eps*u-tau2*par(7)-par(6))/tau1^2;

V_sa= par(1);
P1= par(2);
P2 = par(3);
Cav=par(8);

xq= par(4);
xc= par(5);
xm= par(6);
xm1= par(7);
par=[V_sa;P1;P2;Cav;xq;xc;xm;xm1];

dvdt=[dV_sa(i);dP1(i);dP2(i);dCa(i);dxq(i);dxc(i);dxm1(i);dxm(i)];
dvdt_tot(i)=dvdt;

[Torsten]

I must admit that I don't understand what you are trying to do.

ode45 passes "ABP" and "par" to the function "first".

Here, ABP is a scalar, namely the actual time when your differential equations are to be evaluated.

So what are you trying to do when you generate a new array ABP in "first" and loop over its elements ?

If you try to solve the ODE with a set of different parameters named "ABP", then don't make the loop in "first", but call the integrator ODE45 length(ABP)-times:

ABP= linspace(40,170,131);  
for i=1:length(ABP)
  [T,Y] = ode45(@(t,y)ABP(t,y,ABP(i)),[0 100], [12 97.6 49.67 0.205 0 0 0 0]);
end

function dvdt = first(t,par,abp)
...

Best wishes

Torsten

[Star Strider]

What are you doing in the loop? I do not understand what you want to do.

Your ‘dvdt’ variable should be a single column vector.

Also, this assignment:

dvdt_tot(i)=dvdt;

will fail, because you are assigning a vector to a scalar array element.

[gorilla3]

The aim of the loop is to simply evaluate the equation system for each value of ABP, where ABP is a parameter that ranges from 40 to 170.

[gorilla3]

What you suggested, about moving the loop to the ode45 section makes a lot of sense. Thank you for that!

However, when I remove the outer function (i.e. second), the code does not run and the following error appears:

Subscript indices must either be real positive integers or logicals.

Error in CBF_2006_f_v8>@(t,y)ABP(t,y,ABP(i))

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in CBF_2006_f_v8 (line 42)
  [T,Y] = ode45(@(t,y)ABP(t,y,ABP(i)),[0 100], [12 97.6 49.67 0.205 0 0 0 0]);

[Torsten]

If you save the complete code in one file, the first part also has to be included in a function. So start your file with

function main

ABP= linspace(40,170,131);  
for i=1:length(ABP)
  [T,Y] = ode45(@(t,y)ABP(t,y,ABP(i)),[0 100], [12 97.6 49.67 0.205 0 0 0 0]);
end

function dvdt = first(t,par,abp)
...

and name your file "main.m".

Best wishes

Torsten.

[gorilla3]

Doing that gives me the error: "unable to define function 'main' because it has the same name as the file"