N=200000;
V=zeros(N,1);
V2=zeros(N,1);

G=0.1;
C=1;
E=-70;
Iext=10;
θ=-50
Tarp=1
R=10
Δt=0.01;
t0= t=1 #when I added this command line I obtained multiple spikes. t0 needs to return to the value of t=1.
V[1]=-70;



for k=2:N,

V2[k]= V[k]-50
#V[k]= (Δt/C) * (G*(E-V[k-1]+ Iext)) + V[k-1];

V[k]=V[k-1]+Δt/(R*C) *(E-V[k-1])+(Δt/C)*Iext




        if (V[k] >= θ)
            V[k]=E
            t0=Δt*k
        elseif (Δt*k-t0)<Tarp
            V[k]=E
        else
            V[k]=V[k-1]+Δt/(R*C) *(E-V[k-1])+(Δt/C)*Iext
        end
end


max = (N-1) * Δt;          # Value corresponding to the maximal integration step;
t   = 0:Δt:max;            # Init of an array of N elements, with increasing values (i.e. horizontal coordinate);

using PyPlot;

fig = figure("Action potential",figsize=(12,6));
plot(t, V, "r-", label="membrane potential");
plot(t, V2, "m--", label="membrane potential");
ax = gca();
ax[:set_xlim]((0,50));


ax[:set_ylim]((-100,60));
grid("on");
#xlabel("time (ms)");
#ylabel("mV");
title("Action Potential");
legend(borderaxespad=0.5);