function [Vs,Ws,lams,VRAY,nullV,nullw]=GFitz_Nag_Jac
%  Carol Lucas
%  Gets nullclines, crossing points and eigenvalues (lams) of system

global Af betaf deltaf Cf epsf gamf tend vinitf winitf

clear lams

syms V w

Vdot=Af*V*(V-betaf)*(deltaf-V)-Cf*w;
wdot=epsf*(V-gamf*w);

wsubs=solve(wdot,w);
vsubs=solve(Vdot,w);

[Vs,Ws]=Intersect(Vdot,wdot);

   good=imag(Vs)==0;     %% we only want the real ones (without I, always origin).
   Vs=Vs(find(good));
   Ws=Ws(find(good));

a11=diff(Vdot,V);
a12=diff(Vdot,w);
a21=diff(wdot,V);
a22=diff(wdot,w);
A=[a11 a12;a21 a22];

for k=1:length(Ws)
    w=Ws(k);
    V=Vs(k);
    AN=subs(A);
    lams(:,k)=eig(AN); 
end

VRAY=[-75:.5:75];

for k=1:length(VRAY)
   V=VRAY(k);
 	nullw(k)=subs(wsubs);
	nullV(k)=subs(vsubs);
end

clf;plot(VRAY,subs(nullw),VRAY,subs(nullV),Vs,Ws,'r*');axis([-75 75 -40 40])
xlabel('voltage(MV)');ylabel('w');

LINT=['V = ';'w = '];LEIG=['  eig 1 = ';'  eig 2 = ']; 

for k=1:length(Ws)
    P=num2str([Vs(k);Ws(k)]);
    text(Vs(k),Ws(k),strcat(LINT,P,LEIG,num2str(lams(:,k))));   %% puts info on graph
end
pause(5);