%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global gnabar gkbar Vna Vk Vl Inp BCL 

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

Input = Inp/cosh(mod(t,BCL)-BCL/4);

AM=.1*(25.-V)./(exp(.1*(25.-V))-1.);
BM=4.*exp(-V/18);
AH=.07*exp(-V/20.);
BH=1./(exp(.1*(30.-V))+1.);
AN=.01*(10.-V)./(exp(.1*(10.-V))-1.);
BN=.125*exp(-V/80.);

Fm=AM*(1.-m)-BM*m;
Fh=AH*(1.-h)-BH*h;
Fn=AN*(1.-n)-BN*n;
gna=gnabar*m.^3.*h;
gk=gkbar*n.^4;
Icl=.3*(V-10.5988);
Minf = AM./(AM+BM);
Fv = -gna*(V-Vna)-gk*(V-Vk)-Icl;

s_prime = [Fv+Input Fm Fh Fn]';
