function [Cs,Ws,lams,CRAY,nullC,nullw]=klclosejac(CT)

%% Gets nullclines (nullC, nullw), crossing points(Cs,Ws) and eigenvalues (lams) of system
%% Called by F5_7A

clear lams

syms C w

fi=0.01;kserca=0.2;vserca=100;vRyR=5;vleak=.2;sig=.02;
ka=0.4;kb=0.6;kc=.1;kcm=.1;

winf=(1+(ka/C)^4+(C/kb)^3)/(1+(1/kc)+(ka/C)^4+(C/kb)^3);
CER=(CT-C)/sig;
tau=winf/kcm;
Po=w*(1+(C/kb)^3)/(1+(ka/C)^4 + (C/kb)^3);

jRyR=vRyR*Po*(CER-C);
jserca=vserca*C^2/(kserca^2+C^2);
jleak=vleak*(CER-C);

Cdot=fi*(jRyR+jleak-jserca);
wdot=-(w-winf)/tau;

%% This solves for w in terms of C
wsubs=solve(wdot,w);
csubs=solve(Cdot,w);

%%  Make Cdot a function of only C by substituting for w
Cdote1=subs(Cdot,w,wsubs);

%%  Find all possible values of C
Cs=eval(solve(Cdote1,'C'));

%%  Pick out the real positive ones

good=and(angle(Cs)==0,real(Cs)>=0);
Cs=Cs(find(good));

%%  Now find the Ws that go with those Cs

Ws=subs(wsubs,C,Cs);

a11=diff(Cdot,C);
a12=diff(Cdot,w);
a21=diff(wdot,C);
a22=diff(wdot,w);
A=[a11 a12;a21 a22];

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

CRAY=[0.01:0.01:3];

for k=1:length(CRAY)
   C=CRAY(k);
 	nullw(k)=subs(wsubs);
	nullC(k)=subs(csubs);
end

plot(CRAY,nullw,CRAY,nullC,Cs,Ws,'r*');axis([0 3 0 3])
xlabel('[Ca2+]');ylabel('w');
for k=1:length(Ws)
    text(Cs(k),Ws(k),num2str(lams(:,k)));   %% puts eigenvalues on graph
end
title(['CT = ' num2str(CT)]);
mess='nullcline analysis, hit return to continue',pause;