%     beeler-reuter


Vk=-90;
Vna=50;
V = Vk + (Vna-Vk)*[0:100]/100;
AM=-(V+47)./(exp(-.1*(V+47))-1);
BM=40*exp(-0.056*(V+72));
AH=0.126*exp(-0.25*(V+77));
BH=1.7./(exp(-0.082*(V+22.5))+1);
AJ=0.055*exp(-0.25*(V+78))/(exp(-0.2*(V+78))+1.);
BJ=0.3./(exp(-.1*(V+32))+1.);
Ym=AM./(AM+BM);
Yj = AH./(AH+BH);
Gna = (4*Ym.^3.*Yj.^2+.003);

%Luo-Rudy

V1= Vk +(-40-Vk)*[0:50]/50;
AH1 = 0.135*exp((80+V1)/-6.8);
BH1 = 3.56*exp(0.079*V1)+3.1E5*exp(0.35*V1);
AJ1 = (-1.2714E5*exp(0.2444*V1)-3.474E-5*exp(-0.04391*V1)).*(V1+37.78) ...
  ./(1+exp(0.311*(V1+79.23)));
BJ1 = 0.1212*exp(-0.01052*V1)./(1+exp(-0.1378*(V1+40.14)));

V2 = -40 + (Vna+40)*[1:90]/90;
AH2 = zeros(1,90);
AJ2 = zeros(1,90);
BH2 = 1./(0.13*(1+exp((V2+10.66)/-11.1)));
BJ2 = 0.3*exp(-2.535E-7*V2)./(1+exp(-0.1*(V2+32)));
AH=[AH1,AH2];
BH = [BH1,BH2];
AJ = [AJ1,AJ2];
BJ = [BJ1,BJ2];

V = [V1,V2];
AM=0.32*(V+47.13)./(1-exp(-0.1*(V+47.13)));
BM = 0.08*exp(-V/11);

Zm = AM./(AM+BM);
Zj = AJ./(AJ+BJ);
Zh = AH./(AH+BH);

Tm  =1./(AM+BM);
Tj = 1./(AJ+BJ);
Th = 1./(AH+BH);

plot(V,Tj,'r')




