function [out1, out2, out3] = pNasKModel(t,x,flag)
global E_Na E_K E_Cl V1 V2 i1 i2 R el C Iin eps
V_C=x(1);I_A=x(2);I_S=x(3);I_K=x(4);
if nargin < 3 | isempty(flag)
   out1 = [-(Na(V_C-E_Na)+I_K+Cl(V_C-E_Cl)+I_A-Iin)/C; 
           el*I_S*(V_C-R*I_A);
           el*I_A*(V_C-R*I_A);
           (V_C-K(I_K)-E_K)/eps];
       
else
   switch(flag)
   case 'init'
     out1 = [0; 50];
     out2 = [-0.2507,   -0.3112,    2.7964,   -0.1418];
     out3 = []; %odeset('Vectorized', 1);
     case 'events' 
     %sec1 = V_C-R*I_A; 
     sec1 = V_C-(E_Na+V1);   
     sec2 = V_C-(K(i1)+K(i2))/2-E_K;
     sec3 = V_C-(E_Na+V1); 
     out1 = [sec1,sec2,sec3];
     out2 = [0,0,1];
     out3 = [1,-1,-1];
   otherwise
      error(['Unknown flag ''' flag ''.']);
   end
end; 
