 m    = W(1);
     h    = W(2);
     j    = W(3);
     d    = W(4);
     f    = W(5);
     x    = W(6);
     Cai  = W(7); 
  AM=-(V+47)/(exp(-.1*(V+47))-1);
  BM=40*exp(-.056*(V+72));
  AH=.5*.126*exp(-.25*(V+77));  %should be 77
  BH=.5*1.7/(exp(-.082*(V+22.5))+1);
  AJ=.055*exp(-.25*(V+78))/(exp(-.2*(V+78))+1);%should be 78, not 58
  BJ=.3./(exp(-.1*(V+32))+1);
  AD=.095*exp(-.01*(V-5))./(exp(-.072*(V-5))+1);
  BD=.07*exp(-.017*(V+44))./(exp(.05*(V+44))+1);
  AF=.012*exp(-.008*(V+28))./(exp(.15*(V+28))+1);
  BF=.0065*exp(-.02*(V+30))./(exp(-.2*(V+30))+1);
  AX=.0005*exp(.083*(V+50))./(exp(.057*(V+50))+1);
  BX=.0013*exp(-.06*(V+20))./(exp(-.04*(V+20))+1);
 
  A = [AM, AH, AJ, AD, AF, AX];
  B = [BM, BH, BJ, BD, BF, BX];
  Inf = A./(A+B);                        %    y_inf = a_y/(a_y + b_y);
  m = Inf(1);
  Tau = 1./(A+B);                        
  %    t_y  = 1/(a_y + b_y);
  g_Na = 4*m^3*h*j+0.003;                % Na+ conductance
 
  % The reversal potentials (mv).
  v_Na =  50;                            % sodium equilibrium potential
  % Currents
  INa = g_Na*(V-v_Na);
  IK=1.4*(exp(0.04*(V+85))-1)/(exp(0.08*(V+53))+exp(0.04*(V+53))) ...
  +0.07*(V+23)/(1-exp(-0.04*(V+23)));
  IK = IK + g_katp*(V-Vkatp);
  Ix  = 0.8*x*(exp(0.04*(V+77))-1)/exp(0.04*(V+35));
  Is  = 0.09*d*f*(V+82.3+13.0287*log(Cai));
  if t==t_old | floor(t)==nn+1 
   nn=nn+~(t==t_old);
   Na(nn) = INa;
    K(nn)  = IK;
	X(nn)  = Ix;
	SS(nn)  = Is;
	t_old=t;
	end;
	% Membrane capacitance (micro-Farads/cm^2)
	C_m = 1; % HH
	% Exterior stimulus
	%  I = 3*heaviside(t-100)*heaviside(105-t); 
	%  I = 10*(atan((t-50)*10)/pi+1/2)*(atan((55-t)*10)/pi+1/2); 
	BCL = 200;
	if (mod(t,BCL)>50&mod(t,BCL)<55)
	I = 10;
	else
		I=0.;
	end

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%
	%
	% The derivatives of V, m, h, j, d, f and x with Cai as well.  
	%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	V_prime  = -(INa + IK + Ix + Is - I); 
	%  gate_prime = (Inf - gate)./Tau;   
	gate_prime = A.*(1-gate)-B.*gate;
	 Cai_prime = 0.07*(10^(-7) - Cai) - 10^(-7)*Is;
	 s_prime = [V_prime gate_prime Cai_prime]';
	 
