Δt= 0.001;
N= 200000;
C      = 1.;
#gnamax = 120;
gkmax  = 36;
gl     = 0.3;
Ena    = 50.;
Ek     = -77.;
El     = -54.387;
Iext   = 7;


V= zeros(N,1);
m= zeros(N,1);
n= zeros(N,1);
h= zeros(N,1);
t= zeros(N,1);
τm= zeros(N,1);
τn= zeros(N,1);
τh= zeros(N,1);
αm= zeros(N,1);
βm= zeros(N,1);
αn= zeros(N,1);
βn= zeros(N,1);
αh= zeros(N,1);
βh= zeros(N,1);
minfinity= zeros(N,1);
ninfinity= zeros(N,1);
hinfinity= zeros(N,1);
gnamax= zeros(N:1);

gnamax[1:50000]=0.01
gnamax[50001:100000]=40
gnamax[100001:150000]=80
gnamax[150001:200000]=120

for k=2:N,


V[1]=-64.9964;    #Initial condition;
m[1]=0.0530;      #Initial condition;
n[1]=0.3177;      #Initial condition;
h[1]=0.5960;      #Initial condition;



    αm[k] = 0.1 * (V[k-1]+40.) / (1. - exp(-(V[k-1]+40.)/10.));
    βm[k]  = 4. * exp(-0.0556 * (V[k-1]+65.));
    αn[k] = 0.01 * (V[k-1]+55) / (1. - exp(-(V[k-1]+55.)/10.));
    βn[k] = 0.125 * exp(-(V[k-1]+65.)/80.);
    αh[k] = 0.07 * exp(-0.05*(V[k-1]+65.));
    βh[k] = 1. / (1. + exp(-0.1*(V[k-1]+35.)));

    minfinity[k]  = (αm[k]) / (αm[k] + βm[k]);
    τm[k]  = 1 / (αm[k] + βm[k]);
    m[k]= (Δt*-m[k-1] + Δt*minfinity[k]) / (τm[k]) + m[k-1];

    ninfinity[k]  = (αn[k]) / (αn[k] + βn[k]);
    τn[k]  = 1 / (αn[k] + βn[k]);
    n[k]= (Δt*-n[k-1] + Δt*ninfinity[k]) / (τn[k]) + n[k-1];

    hinfinity[k]  = (αh[k]) / (αh[k] + βh[k]);
    τh[k]  = 1 / (αh[k] + βh[k]);
    h[k]= (Δt*-h[k-1] + Δt*hinfinity[k]) / (τh[k]) + h[k-1];

    V[k]= (Δt/C)*(gnamax[k]*m[k-1]^3*h[k-1]*(Ena - V[k-1]) + gkmax*n[k-1]^4*(Ek - V[k-1]) + gl*(El - V[k-1]) + Iext) + V[k-1];

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;                                 # Tells the computer to “add” a package for (later) generating plots;

fig = figure("Numerical method",figsize=(12,6));    # Create a new figure, with desired size;

#plot(t, V, "r-", label="membrane potential");              # Plot in the current figure, f(x) - the mumerical solution;

ax = gca();                                  # Get a "handle" on the current axes;
ax[:set_xlim]((0,200));                        # For axes "ax" set the horizontal limits;
ax[:set_ylim]((-80,60));                     # For axes "ax" set the vertical limits;
grid("on");                                   # "Grid" on;
xlabel("time");# Label of the horizontal axis
ylabel("mV");                          # Label of the vertical axis
title("Action Potential after pharmacology");        # Label of the title, at the top of the figure
legend(borderaxespad=0.5);                     # Add a legend to the current ax