function connor_stevens

clear 
format long

dt=.05;
tmax=200;
% need integer ratio tmax/dt
nsteps=tmax/dt;
if (floor(nsteps)~=nsteps)
  fprintf(' not an integer number of steps nsteps = %g \n',nsteps)
end

tolnewton=1e-12;
nskip=1;
inewtonmax=20;

Iinj0=8.3;
time1=100;
duration1=30;
Iinj1=10;

time2=160;
duration2=30;
Iinj2=0;

Iinj=Iinj0;

% pre-allocate memory
 maxnumsteps=100000;
 inewton(1:maxnumsteps)=0;
 t(1:maxnumsteps)=dt;
 V(1:maxnumsteps)=0; 
 n(1:maxnumsteps)=0; 
 m(1:maxnumsteps)=0; 
 h(1:maxnumsteps)=0; 
 p(1:maxnumsteps)=0; 
 q(1:maxnumsteps)=0; 

 V(1)=0.0039;
 n(1)=.32;
 m(1)=0.053;
 h(1)=.60;
 p(1)=0;
 q(1)=1;
 
 Y(1,:)=[V(1),n(1),m(1),h(1),p(1),q(1)];
 
 k=1;
 
 while t(k) < tmax
  
  if t(k) >= time1
    Iinj=Iinj1;
  end 
  if t(k) >= time1+duration1
    Iinj=Iinj0;
  end 
  if t(k) >= time2
    Iinj=Iinj2;
  end 
  if t(k) >= time2+duration2
    Iinj=Iinj0;
  end 
   
  [Y(k+1,:),inewton(k)]=BE(Y(k,:),dt,Iinj,tolnewton,inewtonmax,nskip);
  t(k+1)=t(k)+dt;
%  fprintf(' %g %g \n',t(k+1),Iinj);
  
  k=k+1;

  if mod(k,1000) == 0
    fprintf('%g steps t = %g dt =%g \n',k,t(k),dt);
  end
  
 end   


numsteps=k

V=Y(:,1);
Vsignchange=V(1:end-1).*V(2:end);
indsignchange=find(Vsignchange<0);
periods=indsignchange(3:end)-indsignchange(1:end-2);
periods=dt*periods

%figure(3);
%  semilogy(t(1:numsteps),dtj(1:numsteps),'ro');
%  xlabel('Time');
%  ylabel('dt');
   
figure(1);
  subplot(3,1,1)
  plot(t(1:numsteps),Y(1:numsteps,1));
  axis([0 tmax -80 50]);
  subplot(3,1,2)
  legend(' voltage');
  plot(t(1:numsteps),Y(1:numsteps,2),'r',t(1:numsteps),Y(1:numsteps,3),'b',...
       t(1:numsteps),Y(1:numsteps,4),'m');
  legend('n','m','h');
  axis([0 tmax 0 1]);
  subplot(3,1,3)
  plot(t(1:numsteps),Y(1:numsteps,5),'k',t(1:numsteps),Y(1:numsteps,6),'g');  
%  plot(t(1:numsteps),Y(1:numsteps,6),'g'); 
  xlabel('Time')
  axis([0 tmax 0 1]);
  legend('p','q');
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Yp = F(Y,Iinj) 
 
V=Y(1);
n=Y(2);
m=Y(3);
h=Y(4);
p=Y(5);
q=Y(6);

Gna=120;
Ena=55;

Gk=20;
Ek=-72.1;

EA=-75;
GA=47;

Cm=1;
Gl=.3;
El=-17;

Vthresh=-29.7;
alpham=0.38*(V-Vthresh)/(1-exp(-0.1*(V-Vthresh)));
betam=15.2*exp(-0.0556*(V+54.7));
mp=alpham*(1-m)-betam*m;
 
alphah=0.266*exp(-0.05*(V+48));
betah=3.8/(1+exp(-0.1*(V+18)));
hp=alphah*(1-h)-betah*h;

Vthresh=-45.7;
alphan=0.02*(V-Vthresh)/(1-exp(-0.1*(V-Vthresh)));
betan=0.25*exp(-0.0125*(V+55.7));
ninfinity=alphan/(alphan+betan);
taun=1/(alphan+betan);
np=(ninfinity-n)/taun;

alphap=2.5/(1+exp((V+20)*(-0.05)));
betap=0.5/(1+exp((V+0)*(0.08)));
pinfinity=alphap/(alphap+betap);
taup=1/(alphap+betap);
pinfinity=(0.0761*exp(0.0314*(V+94.22))/(1+exp(0.0346*(V+1.17))))^(1/3);
taup=0.3632+1.158/(1+exp(0.0497*(V+55.96)));
pp=(pinfinity-p)/taup;

alphaq=0.2/(1+exp((V+80)*(0.1)));
betaq=0.8/(1+exp((V+50)*(-0.05)));
qinfinity=alphaq/(alphaq+betaq);
tauq=1/(alphaq+betaq);
qinfinity=(1/(1+exp(0.0688*(V+53.3))))^4;
tauq=1.24+2.678/(1+exp(0.0624*(V+50)));
qp=(qinfinity-q)/tauq;

Vp=1/Cm*(Gna*m^3*h*(Ena-V)+Gk*n^4*(Ek-V)+Gl*(El-V)+Iinj+GA*p^3*q*(EA-V));

Yp(1)=Vp;
Yp(2)=np;
Yp(3)=mp;
Yp(4)=hp;
Yp(5)=pp;
Yp(6)=qp;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Y,inewton]=BE(Y,dt,Iinj,tol,inewtonmax,nskip)

  U=Y; 
  inewton=1;
  diffy=1;   
    
  iskip=0;
  while diffy > tol
     
    if mod(iskip,nskip) == 0
      jacobian=dG_BE(dt,U,Y,Iinj);
    end
    iskip=iskip+1;
    gbe=G_BE(dt,U,Y,Iinj)';
    dY=jacobian\gbe;
    U=U-dY';
    diffy=norm(dY);
    inewton=inewton+1;
%    disp([inewton(k),diffy]);
    if inewton > inewtonmax 
     disp(' Newton does not converge');
%     [t(k),u,dy,inewton(k)]
      break
    end
    
  end 
  
  Y=U;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function J=dG_BE(dt,U,Yn,Iinj)

lY=length(Yn);

deltaU=.1e-5;

for j=1:lY
    Up=U;
    Up(j)=Up(j)+deltaU;
    Fp=G_BE(dt,Up,Yn,Iinj);
    J(:,j)=(Fp-G_BE(dt,U,Yn,Iinj))/deltaU;
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


