function dy = chaymodel(t,y)
dy = zeros(3,1);    % a column vector
Vk = -75;   Vi =  100;	Vc =  145;
Vl = -40;   gkv=  1700; gl =  7; gkc=  10.7; gi =  1800;
r  =  0.27;	kc =  3.3/18; 

 
 Ainf = 0.1*(25. + y(1))/(1. - exp(-0.1*y(1) -2.5));
 Binf = 4.*exp(-(y(1) + 50.)/18.);
 Minf = Ainf/(Ainf + Binf);

 Ah   = 0.07*exp(-0.05*y(1) - 2.5);
 Bh   = 1./(1.+exp(-0.1*y(1) - 2.));
 Hinf = Ah/(Ah + Bh);
 
 An = 0.01*(20. + y(1))/(1. - exp(-0.1*y(1)-2.));
 Bn = 0.125*exp( -(y(1)+30.)/80.);
 Ninf = An/( An+Bn );

 Tn = 1./( 230*(An+Bn) );
	 
 dy(1) = gi*Minf^3*Hinf*(Vi-y(1)) + gkv*y(2)^4*(Vk-y(1)) + gkc*y(3)/(y(3)+1.)*(Vk-y(1)) + gl*(Vl-y(1));
 dy(2) = (Ninf - y(2))/Tn;
 dy(3) = r*((Minf^3)*Hinf*(Vc-y(1)) - kc*y(3));