function dy = odebeeler(t,y,options,eNa,eK,gNa,gK)
%y0=[-84,0.5e-7,0,0,1,1,0,1];
%[Time,State]=ode23s('beeler',[0,400],[-84,0.5e-7,0,0,1,1,0,1]);
%plot(Time,State(:,1))
V = y(1);
Ca = y(2);
x = y(3);
m = y(4);
h = y(5);
j = y(6);
d = y(7);
f = y(8);

%gK% = 0.35;
%gNa% = 4.0;

%eK% = -85;
%eNa% = 50;

%gates
c4 = 0; c5 = 0;
ax = (0.0005*exp(0.0830*(V+50)) + c4*(V+c5))/(exp(0.057*(V+50))+1);
bx = (0.0013*exp(-0.060*(V+20)) + c4*(V+c5))/(exp(-0.04*(V+20))+1);
dy(3) = ax*(1-x)-bx*x;
am = (0.0000*exp(0.0000*(V+47)) -  1*(V+47))/(exp(-0.1*(V+47))-1);
bm = (40.000*exp(-0.056*(V+72)) + c4*(V+c5))/(exp(0.00*(V+72))+0);
dy(4) = am*(1-m)-bm*m;
ah = (0.1260*exp(-0.250*(V+77)) + c4*(V+c5))/(exp(0.0*(V+77))+0);
bh = (1.7000*exp(0.0000*(V+22.5))+c4*(V+c5))/(exp(-0.082*(V+22.5))+1);
dy(5) = ah*(1-h)-bh*h;
aj = (0.0550*exp(-0.250*(V+78)) + c4*(V+c5))/(exp(-0.2*(V+78))+1);
bj = (0.3000*exp(0.0000*(V+32)) + c4*(V+c5))/(exp(-0.1*(V+32))+1);
dy(6) = aj*(1-j)-bj*j;
ad = (0.0950*exp(-0.010*(V-5))  + c4*(V+c5))/(exp(-0.072*(V-5))+1);
bd = (0.0700*exp(-0.017*(V+44)) + c4*(V+c5))/(exp(0.05*(V+44))+1);
dy(7) = ad*(1-d)-bd*d;
af = (0.0120*exp(-0.008*(V+28)) + c4*(V+c5))/(exp(0.15*(V+28))+1);
bf = (0.0065*exp(-0.020*(V+30)) + c4*(V+c5))/(exp(-0.2*(V+30))+1);
dy(8) = af*(1-f)-bf*f;

%currents
iK = gK*(4*(exp(0.04*(V-eK))-1)/(exp(0.08*(V+53))+exp(0.04*(V+53)))+0.2*(V+23)/(1-exp(-0.04*(V+23))));
iX = x*0.8*(exp(0.04*(V+77))-1)/exp(0.04*(V+35));
gNaC = 0.003;
iNa = (gNa*m^3*h*j+gNaC)*(V-eNa);
eS = -82.3-13.0287*log(Ca); gS = 0.09;
iS = gS*d*f*(V-eS);

dy(2) = -1e-7*iS+0.07*(1e-7-Ca);


if (t>50.0&t<50.5)
 iEx = -100;
 ist = 0;
else
 iEx = 0;
end

dy(1) = -(1/1)*(iK+iX+iNa+iS+iEx);
%[t y(1)]
dy=dy';