%%  Program for Figure 5.12, Li and Rinzel Model
%%  Carol Lucas
function Css=lrclose_intersect(I,L,Ve)

syms C

fi=0.01;Vi=4;P=26640;
Ki=1.0;Ka=0.4;
Ke=0.2;A=0.5;Kd=0.4;
sigma=0.185;Ct=2;
Ce=(Ct-C)/sigma;
hnull=(Kd./(C+Kd))^3;
cnull=(((Ve*C*C/((Ke*Ke+C*C)*(Ce - C))-L)/P))*((I+Ki)*(C+Ka)/(I*C))^3;
eq=cnull-hnull;
AllCss=solve(eq);
AllCss=subs(AllCss);
good=find(imag(AllCss)==0 & real(AllCss>0));
Css=AllCss(good(length(good)));