%% Figure 6_9
%% Morris_Lecar gap junctions
%% Carol Lucas
%% Calls ode function 'mlcdfun' - seems more stable with ode23

clear tdelray;
clear vmax;
tdelray=[-40:5:-15 -14:14 15:5:40];
v1init=-60.9;w1init=.0149;s2init=0.;s1init=0.;
for k=1:length(tdelray)
    tdel=tdelray(k);
    [t,x]=ode23('mlcdfun',[0 190],[v1init;w1init;s2init;s1init],[],tdel);
    v=x(:,1);w=x(:,2);s1=x(:,4);s2=x(:,3);
    gsyne=1.5;
    vsyne=100;
    isyn2=gsyne*s2.*(vsyne-v);
    isyn1=gsyne*s1.*(vsyne-v);
    figure(1);subplot(211);plot(t,v);axis([0 200 -65 50]);
        title(['tdel = ' int2str(tdel)]),
    subplot(212);plot(t,isyn1,t,isyn2);axis([0 200 -50 200]);
        xlabel('time (ms)');ylabel('Isyn');pause(.5);
    vmax(k)=max(v);
end
figure(2);clf;plot(tdelray,vmax,tdelray,vmax,'ro');
xlabel('time delay');ylabel('maximum voltage of target cell')