%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global gnabar Vna gan Van Vk cm

V = s(1);
m = s(2);
n = s(4);
h = s(3);
      am = -0.1*(V+48.)/(exp(-(V+48.)/15.)-1.);
      bm = 0.12*(V+8.)/(exp((V+8.)/5.)-1.);
      an = -0.0001*(V+50.)/(exp(-(V+50.)/10.)-1.);
      bn = 0.002* exp(-(V+90.)/80.);
      ah = 0.17*exp(-(V+90.)/20.);
      bh = 1./(1.+exp(-(V+42.)/10.));

      gna = gnabar*m.^3*h ;
      gk1 = 1.2*exp(-(V+90.)/50.) + 0.015*exp((V+90.)/60.);
      gk2 = 1.2*n.^4;
      cna = gna*(V-Vna);
      ck = (gk1+gk2)*(V-Vk);
      can = gan*(V-Van);
      Fv = -(cna+ck+can)/cm;
      Fm = am*(1.-m)-bm*m;
      Fh = ah*(1.-h)-bh*h;
      Fn = an*(1.-n)-bn*n;

s_prime = [Fv Fm Fh Fn]';
