% This file  looks at periodically forced solutions of FHN .


clear
time=cputime;



moviename = ['MOVIE.qt']
%MakeQTMovie('start', moviename)


L = 2; % length of the domain
d = .01; % diffusion coefficient
eps = .1;
alf = .05;
gamma = 2.7;
ga = 35;


N=100;                   % number of spatial grid points
M = 120;
h = L/N;
dt = h/4;
Vt=0.015*(1:N)-.5;
Vn = -Vt.*(Vt-1).*(Vt-alf);
Wn = Vt/gamma;

scal = d*dt/h^2;
per = 28;  %gives 2:3 response

per = 400;
T = dt*mod([1:M],per);

Input = 2./cosh(40*(T-dt*per/2)); %monophasic input

%Input = 2./cosh(40*(T-dt*per/2))-4./cosh(40*(T-dt*6.3*per/7)); biphasic
%input
%Input = 0.7*ones(M,1);


X = h*(1:N)';
Inx = 1./cosh(4*X);
V=zeros(N,1); 
W = zeros(N,1);

kstep = 10;
t=0;

% uses Crank Nicolson to solve the diffusion equation

% set up matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);

% The algorithm

for n=2:M
for k = 1:kstep
G = eps*(V-gamma*W);

% update W
W = W + dt*G;
      
F  = ga*(-V.*(V-1).*(V-alf)-W);
Fin = Input(n)*Inx;
V_t =  scal/2 * ([V(2:N);V(N)] - 2*V + [V(1);V(1:N-1)]) + dt*(F + Fin);
V = Am\(V+V_t);
t=t + dt
end


  subplot(2,1,1)
  plot(X,V,'b',X,W,'g',X,Fin,'r')
  axis([0 L -.5 1])

  subplot(2,1,2)
  plot(V,V_t,'b',Vt,Vn,'g--')
 % axis([-.5 1 -.1 .2])
axis([-.5 1 -.1 .3])
  xlabel('{\it v}')
  ylabel('{\it w}')

%MakeQTMovie('addframe','moviename')
%MakeQTMovie('quality', 1.0)

  pause(0.05)

end

%MakeQTMovie('finish','moviename')


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

time_steps = M;

