f = 6.6743*10^(-11);

    m1 = 7.34767309*10^(22);
    m2 = 5.97219*10^(24);
    e = 0.0549006; a = 384748;
    p = a*(1 - e^2);
    c = Sqrt[ f*(m1 + m2)*p];
    eq1 = {x''[t] == - ((
         f*(m1 + m2) *x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^3),
            y''[t] == -((
         f*(m1 + m2) *y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^3),
            z''[t] == -((f*(m1 + m2) *z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^3)};
    rul1 = {Subscript[r, j_] -> Sqrt[
        Subscript[x, j][t]^2 + Subscript[y, j][t]^2 +
         Subscript[z, j][t]^2]};
    initial = {x[0] == 1, x'[0] == 0,
       y[0] == 0, y'[0] == 4,
       z[0] == 0, z'[0] == 3};
sol = NDSolve[eq1, initial, {x[t], y[t], z[t]}, {t, 0, 50}];