%     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);
plot(V,Gna)
