% GFitz_Nag
% Carol Lucas 
% FitzHugh-Nagumo class models
% Try gamf=5 (initial - overdamped), =2 underdamped, =1 oscillatory

global Af betaf deltaf Cf epsf gamf tend vinitf winitf

clear t;

C=input('Hit any number to erase figures: ');
if length(C)~=0; 
	figure(1);clf;figure(2);clf;
	end

tstart=0;

k=1;t=[tstart:.25:tend];

xint=[vinitf;winitf];

[t,x]=ode45('GFNfun',t,xint);

figure(1);subplot(211);plot(t,x(:,1));axis([0 tend -75 75]);
xlabel('time (ms)'),ylabel('voltage (mv)');hold on;
subplot(212);plot(t,x(:,2));
%axis([0 tend 0 1]);
xlabel('time (ms)'),ylabel('W=fraction of open channels');hold on;
pause(5);
mess=('Off to do Phase Plane and Nullclines')

%% Plot phase plane and nullclines

figure(3);clf;   %%  Look here for eigenvalues
[Vs,Ws,lams,VRAY,nullV,nullw]=GFitz_Nag_Jac;


figure(2);plot(x(:,1),x(:,2),x(1,1),x(1,2),'md',VRAY,nullw,':',VRAY,nullV,'--');
ylabel('w'),xlabel('v');axis([-80 80 -40 40]);
title('Phase Plane and Nullclines');hold on;

C=input('Hit any number to get vector field:  ');
if length(C)~=0
 	Vray=[-60:10:60];Nray=[-40:10:40];figure(2);
	vector_field('GFNfun',Vray,Nray,.5);
	end;