function [Vs,Ws,lams,VRAY,nullV,nullw]=GMorris_Lecar_Jac

%% Gets nullclines, crossing points and eigenvalues (lams) of system

global Vk Vl Vca gca gk gl V1 V2 V3 V4 c phi tend Vinit Winit I

clear lams

syms V w

minf=.5*(1+tanh((V-V1)/V2));
winf=.5*(1+tanh((V-V3)/V4));
tauw=1/cosh((V-V3)/(2*V4));
Vdot = ( I - gca*minf*(V-Vca)-gk*w*(V-Vk)-gl*(V-Vl))/c;
wdot=phi*(winf-w)/tauw;

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

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

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,nullw,VRAY,nullV,Vs,Ws,'r*');axis([-75 75 0 1])
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);