Δt= 0.001;  #the formula will work with timesteps of 0,001
N= 200000;    #This will be the size of the array and subsequently the amount of integration steps. You divide (x*N)/Δt = amount of N needed for a certain timeframe.
C      = 1.;    # This is the membrane (specific) capacitance [uF/cm^2]
gnamax = 120;   # This is the membrane max (specific) sodium conductance [mS/cm^2]
gkmax  = 36;    # This is the membrane max (specific) potassium conductance [mS/cm^2]
gl     = 0.3;   # This is the membrane max (specific) leak conductance [mS/cm^2]
Ena    = 50.;     # This is the reversal potential for sodium currents [mV]
Ek     = -77.;    # This is the reversal potential for potassium currents [mV]
El     = -54.387; # This is the reversal potential for leak currents [mV]
Iext   = 7;           #I played with the input current to find the current at which I started getting a train of AP's. There is a threshold at which the current injection causes an action potential to fire. This can be done by trial-and-error adjustment of the external current (Iext). This value lies between 2.2 and 2.25 (nA???) and should correlate with the value at which the frequency goes from 0 to any value higher than 0 in the FI-graph;

V= zeros(N,1); #This command will make the array defined by the number N earlier and fill it with zeros;
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); #because alpha and beta are dependent of "V" and "V" is dependent of time, both of these formulas are vectors and should be put in an array. Later on in the formula of alpham you see that alpha is dependent an V[k-1], so it changes according to changes in V. The same is true for all the alpha's and beta's;
β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);

V[1]=-64.9964;    #This is the resting membrane potential?????? Initial condition; These initial conditions were one cause of the error. I put them before the zeros(N,1) arrays which caused them to be overwritten.
m[1]=0.0530;      #Initial condition;
n[1]=0.3177;      #Initial condition;
h[1]=0.5960;      #Initial condition;

for k = 200:300
Iext = 8
end
for k = 1000:1200
Iext = 8
end

for k=2:N,        #This will start the loop at N = 2. So the first value that will be calculated is K=(2 - 1) => k=1 ;


    #because V is unknown and you have starting point V(1) you can calculate the rest;
    α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]); #This is where the error started. No matching method found;
    τ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*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[k]) + V[k-1]; # This is the hodgkin & huxley model for the generation of action potentials

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);

fig = figure("Numerical method",figsize=(12,6));    # Create a new figure, with desired size;
plot(V, Iext, "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 (ms)");# Label of the horizontal axis
ylabel("mV");                          # Label of the vertical axis
title("Action Potential after paired pulse protocol");        # Label of the title, at the top of the figure
#legend(borderaxespad=0.5);                     # Add a legend to the current ax