% This is a simple ode integrator for the rHH equations
 
function rHH_solver
global gnabar gkbar gl Vna Vk Vl Istim BCL rn0 V

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
%parameters

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

rn0 = 0.8;
%rn0 = 0.9;

Istim = 15;
%Istim = 0;
BCL = 0;

% first find the v nullcline by bisection
vlist = [Vk:0.10001:Vna-0.3];
for j = 1:length(vlist)
    V=vlist(j);
    
% guess a lower bound
a = 0;

% guess an upper bound
b = 1;

%this is the answer - call the function bisect(a,b)

n0(j) = bisect(a,b,@feval);

Gt= gate_de(V);
AN = Gt(3);
BN = Gt(4);
ninf(j) = AN/(AN+BN);

end


tstep = .1;
t_end = 40;
V = 0;
h = 0.3177;

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

s0 = [V;h];  % initial values

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


figure(1)
subplot(2,1,1)
plot(T,V)
xlabel('t')
ylabel('V')
% the phase portrait
subplot(2,1,2)
plot(vlist,n0,'r--',vlist,ninf,'g--',V,h,'linewidth',2)  
xlabel('V')
ylabel('n')
%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global gnabar gkbar gl Vna Vk Vl Istim BCL rn0

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

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

n = W;
Gt= gate_de(V);
AM=Gt(1);
BM = Gt(2);
AN = Gt(3);
BN = Gt(4);
 
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+Istimut Fw]';

 function root = bisect(a,b,feval)

ul = a;
fl = feval(ul);
uu = b;
fu = feval(uu);

% we make the assumption, without checking, that 
% fu*fl<0

% if not, the algoritm fails to find a root.

N = 20;  % number of iterates
% the main bisection algorithm
for j = 1:N
u = (ul+uu)/2;
fc = feval(u);
ftest = (fc*fl>=0);
ul = ftest*u+(1-ftest)*ul;
fl = ftest*fc + (1-ftest)*fl;

uu = (1-ftest)*u+ ftest*uu;
fu = (1-ftest)*fc + ftest*fu;
end
root = u;

function f = feval(u)
    global V
% evaluate the function
% put the function of your choice here
s=[V,u];
t=0;
fhh= deRHS(t,s);
f=fhh(1);

function gate = gate_de(V)
                  
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.);

gate = [AM BM AN BN];
