T = 2000; dt = 0.1;                       % time span and step ( ms )
if length(argv) == 2                      % read from command line
    arg_list = argv();                    % and save parameters in a list
    T  = str2num(arg_list{1});            % run octave script by:
    dt = str2num(arg_list{2});            % > octave izhikevich.m 'T' 'dt'
end;
C = 100; vr = -60; vt = -40; k = 0.7;     % parameters
a = 0.03; b = -2;   c = -50; d = 100;     % neocortical pyramidal neurons
vpeak = 35;                               % spike detected
n = round(T/dt);                          % number of simulation steps
v = vr * ones(1, n); u = 0 * v;           % initial values
I = [zeros(1, 0.1*n), 70*ones(1, 0.9*n)]; % pulse of input DC current
t = 0 * v;                                % time steps

% forward Euler method
for i = 1 : n-1

       % membrane potential
       v(i+1) = v(i) + dt/C * ( k * (v(i) - vr) * (v(i) - vt) - u(i) + I(i) );

       % recovery variable
       u(i+1) = u(i) + dt * a * ( b * (v(i) - vr) - u(i) );

       % save time
       t(i) = i * dt;

       if v(i+1)   >= vpeak               % spike detected !
             v(i)   = vpeak;              % padding spike amplitude
             v(i+1) = c;                  % membrane potential reset
             u(i+1) = u(i) + d;           % recovery variable update
       end;

       % output
       printf ("%6.4e  %8.4e %8.4e\n", t(i), v(i+1), u(i+1));
end;

% generate plot
figure;
plot(t, v, t, u); grid;
xlabel('t (in ms)');
ylabel('V ( in mV) and u (dimensionless)');
legend('Potential', 'Recovery variable');
replot
print ("izhikevich.pdf", "-color", "-dpdf");

