I have no particular problems with this code:
u = 1/2;
r1[x_, y_, u_] := Sqrt[(x + 1 - u)^2 + y^2];
r2[x_, y_, u_] := Sqrt[(x - u)^2 + y^2];
Om[x_, y_, u_] := -1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u];
{xsol, ysol} =
NDSolveValue[{x''[t] - 2 y'[t] == -D[Om[x[t], y[t], u], x[t]],
y''[t] + 2 x'[t] == -D[Om[x[t], y[t], u], y[t]],
x[0] == y[0] == 0, x'[0] == y'[0] == 1/2},
{x, y}, {t, 100}, PrecisionGoal -> 30]
ParametricPlot[{xsol[t], ysol[t]},
{t, 0, 100}]
j[t_] := 1/2*xsol'[t]^2 + 1/2*ysol'[t]^2 +
Om[xsol[t], ysol[t], u];
Plot[j[t], {t, 0, 100}]
