function [xdot,xinit,option] = hh( t, x, flag)
% HH: Hodgkin-Huxley (1952) model of squid axon

global	stim
switch flag
case ''
	% state variables
	v=x(1);		% membrane potential (mV)
	m=x(2);		% sodium current activation
	h=x(3);		% sodium current inactivation
	n=x(4);		% potassium current activation
	% input current (uA/cm^2)
	Ie = interp1(stim(:,1),stim(:,2),t);	% linear interpolation

	% membrane capacitance (uF/cm^2)
	C=1;
	% maximum conductances (uS/cm^2)
	gna = 120;		% sodium
	gk = 36;		% potassium
	gl = 0.3;		% leak
	% reversal potentials (mV)
	Ena = 50;		% sodium
	Ek = -77;		% potassium
	El = -54.4;		% leak
	% activation and inactivation rates
	am = 0.1*(v+40)/(1-exp(-(v+40)/10));
	bm = 4*exp(-(v+65)/18);
	ah = 0.07*exp(-(v+65)/20);
	bh = 1/(1+exp(-(v+35)/10));
	an = 0.01*(v+55)/(1-exp(-(v+55)/10));
	bn = 0.125*exp(-(v+65)/80);
	% xdot
	xdot = [ (Ie - gna*m^3*h*(v-Ena) - gk*n^4*(v-Ek) - gl*(v-El))/C;
		am*(1-m)-bm*m;
		ah*(1-h)-bh*h;
		an*(1-n)-bn*n];
case 'init'			% Initialization
	xdot = stim(1,1):0.1:stim(end,1);	% time span
	xinit = [ -65; 0.1; 0.3; 0.3];	% initial v,m,h,n
	option = odeset('MaxStep',10);	% option
end