function val=kepler(t,x)
    r=sqrt(x(1)^2+x(2)^2);
    val=[x(3);x(4);-x(1)/r^3;-x(2)/r^3];
end

