%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global gnabar Vna
V = s(1);
f = s(2);
h = s(3);
p = s(4);
q = s(5);
d = s(6);

%Remark:  I think there is a mistake in ad

%Remark: q is the slowest variable.
%YNI model:
%'The natural period is about 300'
      tm = V+37.;
      am = (tm>0).*tm./(1.-exp(-tm/10.))+...
	   (tm<0).*tm.*exp(tm/10.)./(exp(tm/10.)-1.)+...
	   (tm==0)*10;
      bm = 40.*exp(-0.056*(V+62.));
      ah = 0.001209*exp(-(V+20.)/6.534);
      bh = 1./(exp(-(V/10.+3.)) + 1.);
      ap = 0.009./(1.+exp(-(V+3.8)/9.71)) + 0.0006;
      bp = .000225*(V+40.)./(exp((V+40.)/13.3) -1.);
      aq = .00034*(V+100.)./(exp((V+100.)/4.4)-1.) +.0000495;
      bq = .0005*(V+40.)./(1-exp(-V/6-6.6667)) + .0000845;
      tm = V+35.;
      ad=(tm>0).*tm./(1.-exp(-tm/2.5))+...
	  (tm<0).*tm.*exp(tm/2.5)./(exp(tm/2.5)-1.)+...
     (tm==0).* 2.5;
      ad = (V>0).*(.01045*ad + .03125*V./(1-exp(-V/4.8)))+...
      (V<0).* 0.01045.*ad + .03125*V.*exp(V/4.8)./(exp(V/4.8)-1.)+...
     (V==0).*( .01045*ad +.03125*4.8);

      bd = .00421*(V-5.)./(exp(V/2.5-2.)-1.);
      af = .000355*(V+20.)./(exp((V+20.)/5.633)-1.);
      bf = .000944*(V+60.)./(1.+exp(-(V+29.5)/4.16));


     Fs = af*(1.-f) - bf*f;
     Fh = ah*(1.-h) - bh*h;
     Fp = ap*(1.-p) - bp*p;
     Fq = aq*(1.-q) - bq*q;
     Fd = ad*(1.-d) - bd*d;
	 minf = am./(am+bm);

      cna = gnabar*minf.^3.*h.*(V-Vna);
      ck = 0.7*p*(exp(0.0277*(V+90.))-1.)/(exp(0.0277*(V+40.)));
      ch = 0.4*q*(V+25.);
      cs = 12.5*(exp(V/15.-2.)-1.)*(0.95*d+0.05)*(0.95*f+0.05);
      cl = 0.8*(1.-exp(-V/20.-3.));
    
      Fv = -(cs + cna + ck + ch + cl);

s_prime = [Fv Fs Fh Fp Fq Fd]';
