function dx = blueSky(t,x, params)

if length(params) == 0
    params = [0.456 0.0357];
end

dx = zeros(3,1);
mu = params(1);
epsilon = params(2);

dx(1) = x(1)*(2+mu-10*(x(1)^2+x(2)^2)) +x(3)^2+x(2)^2+2*x(2);
dx(2) = -x(3)^3-(1+x(2))*(x(3)^2+x(2)^2+2*x(2))-4*x(1) +mu*x(2);
dx(3) = (1+x(2))*x(3)^2+x(1)^2-epsilon;