function [xdot,xinit,option] = cwm( t, x, flag)
% CWM: Connor-Walter-McKown (1977) model of crab 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
	a=x(5);		% A current activation
	b=x(6);		% A current inactivation
	% input current (uA/cm^2)
	Ie = interp1(stim(:,1),stim(:,2),t);	% linear interpretation

	% membrane capacitance (uF/cm^2)
	C=1;
	% maximum conductances (uS/cm^2)
	gna = 120;		% sodium
	gk = 20;	%36;		% potassium
	ga = 47.7;		% A current
	gl = 0.3;		% leak
	% reversal potentials (mV)
	Ena = 55;	%50;		% sodium
	Ek = -72;	%-77;		% potassium
	Ea = -75;
	El = -17;	%-54.4;		% leak
	% activation and inactivation rates
	mshft = -5 -5.3;
	am = 0.1*(v+40+mshft)/(1-exp(-(v+40+mshft)/10));
	bm = 4*exp(-(v+65+mshft)/18);
	hshft = -5 -12;
	ah = 0.07*exp(-(v+65+hshft)/20);
	bh = 1/(1+exp(-(v+35+hshft)/10));
	nshft = -5 -4.3;
	an = 0.01*(v+55+nshft)/(1-exp(-(v+55+nshft)/10));
	bn = 0.125*exp(-(v+65+nshft)/80);
	% A current
	ainf = (0.0761*exp((v+94.22)/31.84)/(1+exp((v+1.17)/28.93)))^(1/3);
	taua = 0.3632+1.158/(1+exp((v+55.96)/20.12));
	binf = 1/(1+exp((v+53.3)/14.54))^4;
	taub = 1.24+2.678/(1+exp((v+50)/16.027));
	% xdot
	xdot = [ (Ie-gna*m^3*h*(v-Ena)-gk*n^4*(v-Ek)-ga*a^3*b*(v-Ea)-gl*(v-El))/C;
		3.8*(am*(1-m)-bm*m);	% 3.8: temperature scaling
		3.8*(ah*(1-h)-bh*h);
		3.8/2*(an*(1-n)-bn*n);
		(-a+ainf)/taua;
		(-b+binf)/taub];
case 'init'			% Initialization
	xdot = [stim(1,1) stim(end,1)];	% time span
	xinit = [ -65; 0.1; 0.1; 0.5; 0.5; 0.5];	% initial v,m,h,n,a,b
	option = odeset('MaxStep',1,'RelTol',1e-5);	% options
end