Δt= 0.01;  #the formula will work with timesteps of 0,01
N= 50000;    #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]




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

Iext=7;           #I played with the input current to find the current at which I started getting a train of AP's;

for k=2:N,        #This will start the loop #Not sure what to put at the question mark k=???:N;



    αm[k] = 0.1 * (V[k-1]+40.) / (1. - exp(-(V[k-1]+40.)/10.));  #because V is unknown and you have starting point V(1) you can calculate the rest;
    β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) + V[k-1]; # This is the hodgkin & huxley model

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