function fty = f(t,y)
        fty(1,1) = y(3);
        fty(2,1) = y(4);
        fty(3,1) = -y(1)/(y(1)^2 + y(2)^2)^(3/2);
        fty(4,1) = -y(2)/(y(1)^2 + y(2)^2)^(3/2);
        return;

