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

clear all
global gam alf eps Inp BCL
alf = 0.139;
eps = 0.005;
gam = 2.54;

Inp = .1;
BCL = 155;
tstep = 1;
t_end = 2000;
V = 0;
W = 0;

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

s0 = [V;W];

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

FHN_nullclines

figure(1)
subplot(2,1,1)
plot(T,V,'r',T,W,'b')
subplot(2,1,2)
plot(V,W,Vt,Vn,'b--',Vt,Wn,'r--')
axis(wrange)
xlabel('V')
ylabel('W')

%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global gam alf eps Inp BCL
V = s(1);
W = s(2);

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

Fw = eps*(V - gam*W);
Fv = -V.*(V-1).*(V-alf) - W ;
s_prime = [Fv+Input Fw]';
