% This is a simple ode integrator for the rHH equations
% it is set up to be easily modified to use for other dynamics
function rHH
clear all

global gnabar gkbar gl Vna Vk Vl Inp BCL rn0
gnabar = 120.;
gkbar = 36;
gl = 0.3;
Vna = 115;
Vk = -12;
Vl = 10.5988;

rn0 = 0.85;
%rn0 = 0.9138;
Inp = 5;
BCL = 10;
BCL = 0; % this makes the input constant
tstep = .1;
t_end = 40;
V = 0;
h = 0.3177;

%specify the output points
tspan = [0:tstep:t_end];

s0 = [V;h];

[T,S] = ode23(@deRHS,tspan, s0, odeset('maxstep',1));  
V = S(:,1);
h = S(:,2);

figure(1)
subplot(2,1,1)
plot(T,V)

rHH_nullclines
subplot(2,1,2)
plot(V,h,'r',Vn,n0,'b--',Vn,ninf,'g--')  

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

V = s(1);
W = s(2);

Input = (BCL>0)*Inp/cosh(mod(t,BCL)-BCL/4) + (BCL==0)*Inp;

n = W;
     
AM=.1*(25.-V)./(exp(2.5 -0.1*V)-1.);
BM=4.*exp(-V/18.);
AN=.01*(10.-V)./(exp(1.-0.1*V)-1.);
BN=.125*exp(-V/80.);

h = rn0 - n;
aminf = AM/(AM+BM);
aminf = aminf.^3;
Ina = gnabar*aminf.*h.*(V-Vna) ;
Ik = gkbar*n.^4*(V-Vk);
Icl = gl*(V-Vl);
      
Fv = -Ina-Ik-Icl;
Fw=AN*(1.-n)-BN*n;



s_prime = [Fv+Input Fw]';
