
function [t, V] = HodgkinHuxley()

V0 = -75;
n0 = (0.01*(V0+10)/exp(((V0+10)/10)-1))/((0.01*(V0+10)/exp(((V0+10)/10)-1))+(0.125*exp(V0/80)));
m0 = (0.1*(V0+25)*(exp(((V0+25)/10)-1)))/((0.1*(V0+25)*((exp(((V0+25)/10)-1))))+(4*exp(V0/18)));
h0 = (0.07*exp(V0/20))/((0.07*exp(V0/20))+(1/(exp(((V0+30)/10)+1))));

[t, V] = ode45(@deriv,[0,100],[V0;n0;m0;h0]);

plot(t,V(:,1));



function dxy = deriv(t,xy)

V = xy(1);
n = xy(2);
m = xy(3);
h = xy(4);

Cm = 1.0; % microF per cm^2
VNa = -115; % mV
VK = 12; % mV
VL = -10.613 ; %mV
GNa = 120; % mS per cm^2 (called mmho in paper)
GK = 36; %mS per cm^2
GL = 0.3; %mS per cm^2

if t==0.01
    Istim = 20;
else
    Istim = 0;
end

an = 0.01*(V+10)/exp(((V+10)/10)-1);
bn = 0.125*exp(V/80);
am = 0.1*(V+25)*(exp(((V+25)/10)-1));
bm = 4*exp(V/18);
ah = 0.07*exp(V/20);
bh = 1/(exp(((V+30)/10)+1));


dVdt = (Istim - (GK*(n^4)*(V-VK)) - (GNa*(m^3)*h*(V-VNa)) - (GL*(V-VL)))/Cm
dndt = an*(1-n) - bn*n;
dmdt = am*(1-m) - bm*m;
dhdt = ah*(1-h) - bh*h;

dxy = [dVdt;dndt;dmdt;dhdt];