GK = 36; GNa = 120; VL = 10.6;
 alpham[V_] := (0.1 (25 - V))/(Exp[(25 - V)/10] - 1);
 betam[V_] := 4 Exp[-V/18];
 alphan[V_] := (0.01 (10 - V))/(Exp[(10 - V)/10] - 1);
 betan[V_] := 0.125 Exp[-V/80];
 alphah[V_] := 0.07 Exp[-V/20];
 betah[V_] := 1/(Exp[(30 - V)/10] + 1);
 gNa[t_] := ((m[t])^3) h[t];
 gK[t_] := (n[t])^4;
Istim[t_, duration_, strength_] := If[0 < t < 100, strength, 0];


Manipulate[
Show[
  Plot[
   Evaluate[
    V[t] /. NDSolve[{V'[t] ==
        1/Cm (Istim[t, duration,
            strength] - ((GNa) gNa[t] (V[t] - VNa)) - ((GK) gK[
              t] (V[t] - VK)) - ((GL) (V[t] - VL))),
       m'[t] == (alpham[V[t]] (1 - m[t])) - (betam[V[t]] m[t]),
       n'[t] == (alphan[V[t]] (1 - n[t])) - (betan[V[t]] n[t]),
       h'[t] == (alphah[V[t]] (1 - h[t])) - (betah[V[t]] h[t]),
       V[0] == -0.0741, m[0] == 0.0525, n[0] == 0.3183,
       h[0] == 0.588}, {V, m, n, h}, {t, 0, 80},
      MaxSteps -> 100000]],
   {t, 0, 60}, PlotRange -> {-20, (VNa + 10)}, AspectRatio -> 1,
   PlotLabel -> "membrane voltage versus time",
   AxesLabel -> {"time (msec)", "membrane voltage (mV)"},
   ImageSize -> {550, 300}, PlotStyle -> {Thickness[.015], Blue}],
  Graphics[{Thickness[0.008], Dashed, Red,
    Line[{{0, VNa}, {60, VNa}}], Line[{{0, VK}, {60, VK}}]}]],
{{strength, 6.3, "applied current"}, 0, 18,
  Appearance -> "Labeled"},
{{VNa, 115, "sodium equilibrium potential"}, 110, 150,
  Appearance -> "Labeled"},
{{VK, -12, "potassium equilibrium potential"}, -15, 5,
  Appearance -> "Labeled"},
Delimiter,
{{Cm, 1, "membrane capacitance"}, {1, 10^-2.99}, Checkbox},
{{GL, 0.3, "leakage channel conductance"}, {0.3, 0}, Checkbox},
SaveDefinitions -> True]

 (*Values for the physical quantities relevant in the experiment performed by Hodgkin and Huxley with the squid giant axon. \
 Conductances are given in mho/cm^2, voltages in mV*)
 GNa = 120.0*10^-3; (*Sodium Maximum Conductance*)
 GK = 36*10^-3; (*Potassium Maximum Conductance*)
 GL = 0.3*10^-3; (*Leakage Current Maximum Conductance*)
 VNa = 115; (*Sodium Diffusion Potential*)
 VK = -12; (*Potassium Diffusion Potential*)
 VL = 11; (*Leakage Current Diffusion Potential*)
 a = 238*10^-4; (*Radius of the axon, given in cm*)
Cap = 10^-3; (*Capacitance per unit area of the axon membrane, given in mF/cm^2*)
r = 2*10^4; (*Longitudinal resistance per unit length of the axon, given in ohm/cm*)
m0 = 0.053;
h0 = 0.59;
n0 = 0.318;
V0 = 0;

(* Define utility functions *)

(*Voltage-dependent coefficients to be used in the equations for the dynamics of gate-variables m,h and n*)
am[V_] := (0.1 (25 - V))/(E^((25 - V)/10) - 1);
bm[V_] := 4 E^(-(V/18));
ah[V_] := 0.07 E^(-(V/20));
bh[V_] := 1/(E^((30 - V)/10) + 1);
an[V_] := (0.01 (10 - V))/(E^((10 - V)/10) - 1);
bn[V_] := 0.125 E^(-(V/80));
(*Equations for the dynamics of gate-variables m,h and n. dm, dh and dn stand for dm/dt, dh/dt and dn/dt, respectively*)
dm[t_] := am[V[t]] (1 - m[t]) - bm[V[t]] m[t];
dh[t_] := ah[V[t]] (1 - h[t]) - bh[V[t]] h[t];
dn[t_] := an[V[t]] (1 - n[t]) - bn[V[t]] n[t];
(*Aperture of the sodium and potassium channels*)
ChannelsNa[u_] := (m[u])^3 h[u];
ChannelsK[u_] := (n[u])^4;
(*Current density of sodium, potassium and leakage current, together with the resultant ion current density*)
JNa[u_] := GNa ChannelsNa[u] (V[u] - VNa);
JK[u_] := GK ChannelsK[u] (V[u] - VK);
JL[u_] := GL (V[u] - VL);
Jion[u_] := JNa[u] + JK[u] + JL[u];
(*UPS stands for 'Unique Point Solution', since it finds and plots a numerical solution for the behaviour of voltage, currents and gates in a chosen point along the axon. J0 is the stimulus current,dV is the function dV/dt and s is the numerical solution returned by the Mathematica function NDSolve*)

(* Initialize stimulus and plot range values *)

stimulus = 0.02;(* stimulus strength (mA/cm^2) *)
start = 10;(* start of stimulus (ms) *)
duration = 40;(* stimulus duration (ms) *)
domaini = 0;(* initial time *)
domainf = 50;(* final time *)
rangel = -20;(* minimum range *)
rangeh = 120;(* maximum range *)

(* Numerically solve the DE and plot the potential versus time *)

Block[{J0, dV, s},
J0 = Function[t, If[start <= t <= duration, stimulus, 0]];
dV = Function[t, (J0[t] - Jion[t])/Cap];
s = NDSolve[{
    V'[t] == dV[t],
    m'[t] == dm[t],
    h'[t] == dh[t],
    n'[t] == dn[t],
    V[0] == V0,
    m[0] == m0,
    h[0] == h0,
    n[0] == n0},
   {V, m, h, n},
   {t, domainf}
   ];
Plot[
  Evaluate[{V[t]} /. s],
  {t, domaini, domainf},
  PlotRange -> {rangel, rangeh},
  AxesLabel -> {"time (ms)", "membrane voltage (mV)"},
  AxesOrigin -> {domaini, rangel},
  GridLines -> Automatic, GridLinesStyle -> Dashed,
  ImageSize -> 450, PlotStyle -> {Blue}]
]