% This file  looks at two coupled  FHN cells 
% as a model of border zone arrhythmia.


clear
time=cputime;

eps = .11;
alf = .05;
gamma = 2.5;
beta = .25;
%oscillations occur with coup = 2, but not with coup = 0;
coup = 2;
%coup = 0;

M = 200;
kstep = 20;
dt = 0.01;
N=100;


%phase plane
Vt=0.015*(1:N)-.5;
V1n = -Vt.*(Vt-1).*(Vt-alf);
V2n =-(1-beta)*Vt.*(Vt-1).*(Vt-alf)+beta*(1-Vt);
Wn = Vt/gamma;
V1=0; 
W1 = 0;
v1(1,1)=V1;
w1(1,1)=W1;
V2=0; 
W2 = 0;
v2(1,1)=V2;
w2(1,1)=W2;
Tt(1,1)=0;



t=0;

% The algorithm

for n=2:M
for k = 1:kstep
G1 = eps*(V1-gamma*W1);
G2 = eps*(V2-gamma*W2);

% update W
W1 = W1 + dt*G1;
W2 = W2 + dt*G2;
      
F1  = 25*(-V1.*(V1-1).*(V1-alf)-W1) +coup*(V2-V1);
      
F2  = 25*(-(1-beta)*V2.*(V2-1).*(V2-alf)+beta*(1-V2)-W2) +coup*(V1-V2);

V1 = V1+dt*F1;
V2=V2+dt*F2;

t=t + dt;
v1(n,1)=V1;
w1(n,1)=W1;
v2(n,1)=V2;
w2(n,1)=W2;
Tt(n,1)=t;
end
end

%for n=1:M
figure(1);
title(['coupled cells']);

  subplot(2,2,1)
  plot(Tt,v1,'b',Tt,w1,'g',Tt(n,1),v1(n,1),'*',Tt(n,1),w1(n,1),'*')
  title(['normal cell']);
  subplot(2,2,2)
  plot(Tt,v2,'b',Tt,w2,'g',Tt(n,1),v2(n,1),'*',Tt(n,1),w2(n,1),'*')
  title(['ischemic cell']);
  
  subplot(2,2,3)
  plot(v1(n,1),w1(n,1),'*',Vt,V1n,'g--',Vt,Wn,'r--',v1,w1,'b')
  axis([-.5 1 -.1 .3])
  
   subplot(2,2,4)
  plot(v2(n,1),w2(n,1),'*',Vt,V2n,'g--',Vt,Wn,'r--',v2,w2,'b')
  axis([-.5 1 -.1 .3])

  pause(0.05)

%end

TL = t;
time=cputime-time        % outputs the cpu time taken to solve the equations

time_steps = M;

