How to define initial conditions without getting an array error

I need to define some initial conditions but I keep getting an "Array indices must be positive integers or logical values" error. I don't really know how to fix this... also the initial conditions must be defined in a column vector. I've written the initial conditions that I was given as a side note in my code. Thank you!

Here is my code:

clear;

clc;

timespan=[0 170]; % -20 from stim on (20-20); +20 from stim off (150 + 20)...

c=struct;

c.cm=1;

c.vrest=-70;

c.gca=1.1;

c.gk=2;

c.eca=100;

c.ek=-70;

c.gm=0.5;

vm=-30; % initial condition: Vm(0)=-30

w=0.01; % initial condition: w(0)=0.1

y_0=[vm(0);w(0)]; % has to be in column vector!

options=[];

[t,y]=ode45(@odefun,tspan,y_0,options,c);

function [output]=odefun(t,y,c)

w_infnty=(1+tanh(y(1)/30))/2;

m_infnty=(1+tanh(y(1)+1)/15)/2;

T_w=5/cosh(y(1)/60);

if t>=20 && t<= 150

pulse_t=40;

else

pulse_t=0;

end

dy_dt(1)=-c.gm*(y(1)-c.vrest)/c.cm-c.gca*m_infnty*(y(1)-c.eca)/c.cm-c.gk*y(2)*(y(1)-c.ek)/c.cm+pulse_t/c.cm;

dy_dt(2)=(w_infnty-y(2))/T_w;

output=[dy_dt(1); dy_dt(2)];

end

Answer 1

[Sam Chak]

Hi @Gabriela,

The initial values can be specified as scalar values, without using array indexing method.

tspan   = [0 30];   % -20 from stim on (20-20); +20 from stim off (150 + 20)...
c       = struct;
c.cm    = 1;
c.vrest = -70;
c.gca   = 1.1;
c.gk    = 2;
c.eca   = 100;
c.ek    = -70;
c.gm    = 0.5;
vm0     = -30;          % initial condition: Vm(0)=-30
w0      = 0.01;         % initial condition: w(0)=0.1
y_0     = [vm0; w0];    % has to be in column vector!
% options = [];
[t, y]  = ode45(@(t, y) odefun(t, y, c), tspan, y_0);
plot(t, y), grid on
xlabel('t'), ylabel('\bf{y}')
legend('y_1', 'y_2')

function output = odefun(t, y, c)
    w_infnty = (1 + tanh(y(1)/30))/2;
    m_infnty = (1 + tanh(y(1) + 1)/15)/2;
    T_w      = 5/cosh(y(1)/60); 
    if t>=20 && t<= 150
        pulse_t = 40;
    else
        pulse_t = 0;
    end 
    dy_dt(1) = -c.gm*(y(1)-c.vrest)/c.cm-c.gca*m_infnty*(y(1)-c.eca)/c.cm-c.gk*y(2)*(y(1)-c.ek)/c.cm+pulse_t/c.cm;
    dy_dt(2) = (w_infnty-y(2))/T_w;
    output   = [dy_dt(1); dy_dt(2)];
end