% GMorris_Lecar_Bif.m
% Carol Lucas
% Simulations for Figure 2.9D
% Parameters come from Table 2.4
% Parameters should be initialized with GMorris_Lecar_Init.m
% This is brute force - takes a long time

global Vk Vl Vca gca gk gl V1 V2 V3 V4 c phi tend Vinit Winit 
syms V w I
tstart=0;tend=800;Isave=I;

k=1;t=[tstart:2:tend];

clear Vmin Vmax Vest;

xint=[Vinit;Winit];
k=1;IRAY=[0:10:90 91:99 100:10:230 231:239 240:10:300];

minf=.5*(1+tanh((V-V1)/V2));
winf=.5*(1+tanh((V-V3)/V4));
tauw=1/cosh((V-V3)/(2*V4));
wdot=phi*(winf-w)/tauw;

tic
for k=1:length(IRAY);
	I=IRAY(k);   
	[t,x]=ode45('Gmlfun',t,xint);
	Vss=x([length(t)-100:length(t)],1);
	Vmin(k)=min(Vss);
	Vmax(k)=max(Vss);
    Vdot = ( I - gca*minf*(V-Vca)-gk*w*(V-Vk)-gl*(V-Vl))/c;
   [Vs,Ws]=Intersect(Vdot,wdot);
   Vest(k)=Vs;
   West(k)=Ws;
    P=strcat('I = ',num2str(I),'  Split between min and max = ',num2str(Vmax(k)-Vmin(k),3))
  

end;
toc
figure(3);clf;plot(IRAY,Vmin,IRAY,Vmax,IRAY,Vest,'r:');axis([0 300 -75 75]);
xlabel('Iapp'),ylabel('voltage (mv)');
title('BIFURCATION ANALYSIS');
I=Isave;  %% Reset to original I