function dx = hr( t,x,params )

if length(params) == 0
    params = [2.96 2.99 0.01 4];
end

b = params(1);
I = params(2);
mu = params(3);
s = params(4);

dx = zeros(3,1);

dx(1) = x(2) - x(1)^3 + b*x(1)^2 + I - x(3);
dx(2) = 1 - 5*x(1)^2 - x(2);
dx(3) = mu*(s*(x(1) + 1.6) - x(3));
%mu = 0.002;
%dx(3) = mu*(s*(x(1)+1.4)-x(3)-0.0083/((x(3)-0.9)^2+0.03));

end

