function dydt = f(t,y,flag, Xbar, k)
%Xbar = [15000 15000 15000]/2^16 ;
%k = [656 52224 32768 10496 8192]/2^16 ;
% [S] = y(1), [ES] = y(2), [P] = y(3)
dydt = [k(1)*Xbar(1)*y(2) - k(2)*y(1)*y(2) + k(3)*Xbar(2)*y(1) - 2*k(4)*y(1)*y(1) ;...
       -k(1)*Xbar(1)*y(2) - k(2)*y(1)*y(2) + k(5)*Xbar(3)*y(3) ; ...
        k(3)*Xbar(2)*y(1) - k(5)*Xbar(3)*y(3) ; ...
        ];