function [x y u w] = Apollo(x0, y0, u0, w0, f, t, n)
    h = t / n;
    x = [ x0 ]; y = [ y0 ]; u = [ u0 ]; w = [ w0 ];

    m1 = 1 / 81.45;
    m2 = 1 - m1;

    r1 = @(x, y) ((x + m1)^2 + y^2) ^ (3/2);
    r2 = @(x, y) ((x - m2)^ 2 + y^2) ^ (3/2);

    f1 = @(x, y, u, w) 2*w + x - f*u - m2*(x+m1)/r1(x,y) - m1*(x-m2)/r2(x,y);
    f2 = @(x, y, u, w) -2*u + y - f*w - m2*y/r1(x,y) - m1*y/r2(x,y);

    for i = 1 : n
        kx1 = h * u(i);
        ky1 = h * w(i);
        ku1 = h * f1(x(i), y(i), u(i), w(i));
        kw1 = h * f2(x(i), y(i), u(i), w(i));

        kx2 = h * (u(i) + ku1/2);
        ky2 = h * (w(i) + kw1/2);
        ku2 = h * f1(x(i)+kx1/2, y(i)+ky1/2, u(i)+ku1/2, w(i)+kw1/2);
        kw2 = h * f2(x(i)+kx1/2, y(i)+ky1/2, u(i)+ku1/2, w(i)+kw1/2);

        kx3 = h * (u(i) + ku2/2);
        ky3 = h * (w(i) + kw2/2);
        ku3 = h * f1(x(i)+kx2/2, y(i)+ky2/2, u(i)+ku2/2, w(i)+kw2/2);
        kw3 = h * f2(x(i)+kx2/2, y(i)+ky2/2, u(i)+ku2/2, w(i)+kw2/2);

        kx4 = h * (u(i) + ku3);
        ky4 = h * (w(i) + kw3);
        ku4 = h * f1(x(i)+kx3, y(i)+ky3, u(i)+ku3, w(i)+kw3);
        kw4 = h * f2(x(i)+kx3, y(i)+ky3, u(i)+ku3, w(i)+kw3);

        x(i+1) = x(i) + (kx1 + 2*kx2 + 2*kx3 + kx4) / 6;
        y(i+1) = y(i) + (ky1 + 2*ky2 + 2*ky3 + ky4) / 6;
        u(i+1) = u(i) + (ku1 + 2*ku2 + 2*ku3 + ku4) / 6;
        w(i+1) = w(i) + (kw1 + 2*kw2 + 2*kw3 + kw4) / 6;
    end
end