function dy = FitzhughNagumo_eqn(t,y)
global a epsi gamma
dy = zeros(2,1); % need to be a column vector
dy(1) = -y(1)*(y(1)-a)*(y(1)-1)-y(2);
dy(2) = epsi*(y(1)-gamma*y(2));