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


%parameter values for the Mitchell-Schaeffer model
global topen tclose gna gk vh vna vk kh Inp BCL
tclose = 150;
topen = 120;
%tclose = 200;
gna = 3.33; 
%gna = 2.;
gk = 0.167;
gk = .1;
vh = 0.13;

vna = 1;
vk = 0.012;
kh = .02; 

%cardiac_parameters %take your pick of parameter values
Inp = .1;
Inp = 0;
BCL = 372;
tstep = 2;
t_end = 2000;
V = 0;
W = .2;

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

s0 = [V;W];

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

MS_nullclines

tplot = [1:t_end/tstep];
figure(1)
subplot(2,1,1)
plot(T(tplot),V(tplot),'r',T(tplot),W(tplot),'b','linewidth',2)
xlabel('time','fontsize',16)
ylabel('\phi, {\it h}','fontsize',16)
subplot(2,1,2)
plot(V(tplot),W(tplot),'linewidth',2)
hold on
plot(Vt,Vn,'b--',Vt,Wn,'r--')
hold off
axis(wrange)
xlabel('\phi','fontsize',16)
ylabel('{\it h}','fontsize',16)


%the right hand side for ode simulation:
function s_prime=deRHS(t,s)
global topen tclose gna gk vh vna vk kh Inp BCL

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

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


minf = max(0,V);
finf = 0.5*(1.+tanh((V-vh)/kh));
   hinf = 1-finf;
   tauh = topen+(tclose-topen)*finf;
   Fv = gna*W.*max(0,V).^2.*(vna-V) +gk*(vk-V);
   Fw = (hinf-W)./tauh;

s_prime = [Fv+Input Fw]';

