Error using odearguments (line 93); must return a column vector.

Hi guys!

I'm trying to run a stimulation of two neural masses in matlab, however the ode is giving me a error code.

I tried to fix it myself but after some trial and error I keep getting the same error.

Down below is both my run script, which runs the ode, and the function itself.

%Firing rate equations for 2 neural masses with only electrical coupling
%& chemical coupling. N=10000 

I_1= pi^2;
I_2= -2*pi^2;

params.tau= 10;
params.g= 1;
params.delta= 1;
N= 10000;


V_1= ones(10000,1)*-65;
R_1= zeros(10000,1);
V_2= ones(10000,1)*-65;
R_2= zeros(10000,1);
I= [I_1;I_2];

y0=[V_1,V_2,R_1,R_2];

options = odeset('Events',@QIF_reset)


tstart=0
tfinal= 100
tout=tstart
yout = y0.';
teout = [];
yeout = [];
ieout = [];


for i = 1:5
   % Solve until the first terminal event.
   [t,y,te,ye,ie] = ode23(@(t,y)QIF(t,y0,I,params,N),[tstart tfinal],y0,options);
   if ~ishold
      hold on
   end
   % Accumulate output.  This could be passed out as output arguments.
   nt = length(t);
   tout = [tout; t(2:nt)];
   yout = [yout; y(2:nt,:)];
   teout = [teout; te];          % Events at tstart are never reported.
   yeout = [yeout; ye];
   ieout = [ieout; ie];
   

   % Set the new initial conditions, with .9 attenuation.
   y0(:,1) = -100;
   y0(:,2) = -65;
   y0(:,3) = 0;
   y0(:,4) = 0;
   
   % A good guess of a valid first timestep is the length of the last valid
   % timestep, so use it for faster computation.  'refine' is 4 by default.
  
   tstart = t(nt);
end
plot(tout,yout(:,1))

And here is the function which the run script runs:

function yp= QIF(~,y,I,params,N)
%QIF ordinary differential equation solver for two neurons electrically
%coupled. The input of these neurons are n1= pi^2 and n2= -2*pi^2
I_1= I(1);
I_2=I(2);
v_1= sum(y(1:10000,1))/N;
v_2= sum(y(1:10000,2))/N;

for i=1:N
    
V_1d(i)= (((y(i,1).^2)+ I_1)+params.g*(v_1-y(i,1)))/params.tau;
V_2d(i)= (((y(i,2).^2)+ I_2)+params.g*(v_2-y(i,2)))/params.tau;;

R_1d(i)= ((params.delta/(params.tau*pi))+ 2*y(i,3)*y(1)-params.g*y(i,3));
R_2d(i)= ((params.delta/(params.tau*pi))+ 2*y(i,4)*y(2)-params.g*y(i,4));
end

yp = [V_1d;V_2d;R_1d;R_2d];

Answer 1

[Steven Lord]

Well, does your function return a column vector like the error message states it must, or does it return a row vector?

The easiest solution will probably be to reshape yp to be a column vector before returning it from your function.

[Maarten Duijn]

Yes, that's the part I understand obviously.

However, I tried both yp= [x;y;z] and yp= [x,y,z] both give the same error code.

How would you suggest reshaping yp otherwise?

[Steven Lord]

yp = [V_1d;V_2d;R_1d;R_2d];
yp = reshape(yp, 4*N, 1);

Be careful about how your initial conditions are ordered, to ensure they're ordered in the same order as this yp vector. Remember that MATLAB will reshape the matrix in column major order.

V_1d = 1:4;
V_2d = 5:8;
R_1d = 9:12;
R_2d = 13:16;

yp = [V_1d;V_2d;R_1d;R_2d];
yp = reshape(yp, 4*N, 1) % This is not 1:16

If you expected yp to be 1:16, you probably want to preallocate each of your vectors to be column vectors before you fill them in.

V_1d = zeros(N, 1); % and the same for V_2d, R_1d, R_2d

If you preallocate them to be column vectors, then you won't need to use reshape.

V_1d = zeros(N, 1);
V_1d(1) = 1;
V_1d(2) = 4;
V_1d(3) = 9;
V_1d(4) = 16;

V_2d = V_1d + 1;
R_1d = V_1d + 2;
R_2d = V_1d + 3;

yp = [V_1d;V_2d;R_1d;R_2d]

Answer 2

[Marcel Jiokeng]

You should initialize yp in the function QIF(~,y,I,params,N) as it follows:

function yp = QIF(~,y,I,params,N)

yp = zeros((with the right length),1);

....

....

....