%this file finds the nullclines for the rHH equations

Vn =.01+ [-10:20]';

AM=.1*(25.-Vn)./(exp(2.5 -0.1*Vn)-1.);
BM=4.*exp(-Vn/18.);
AN=.01*(10.-Vn)./(exp(1.-0.1*Vn)-1.);
BN=.125*exp(-Vn/80.);

ninf = AN./(AN+BN);
aminf = AM./(AM+BM);
aminf = aminf.^3;

%make this a two component problem
Nv = length(Vn);
na = zeros(Nv,2);
n0 = na;
eval_Fvn
Fa = Fvn;


nb = ones(Nv,2);
n0 = nb;
eval_Fvn
Fb = Fvn;
%plot(Vn,Fa(:,1))
%hold on
%plot(Vn,Fb(:,1))

for j = 1:15
    nc = (na+nb)/2;
    n0 = nc;
    eval_Fvn
    Fc = Fvn;
    tmp = Fc.*Fb;
    [indx,jndx] = find(tmp>0);
   Fb(indx,jndx) = Fc(indx,jndx);
    nb(indx,jndx) = nc(indx,jndx);
    [kndx,lndx] = find(tmp<=0);
    Fa(kndx,lndx) = Fc(kndx,lndx);
    na(kndx,lndx) = nc(kndx,lndx);
  
hold on

end

hold off

