def METHOD_NAME(Y):
    # defining the ODE:
    # assume Y = [y1,y2,y1',y2']
    mu = 0.012277471  # a constant (standardised moon mass)
    mp = 1 - mu
    D1 = ((Y[0] + mu) ** 2 + Y[1] ** 2) ** (3 / 2)
    D2 = ((Y[0] - mp) ** 2 + Y[1] ** 2) ** (3 / 2)
    y1p = Y[0] + 2 * Y[3] - mp * (Y[0] + mu) / D1 - mu * (Y[0] - mp) / D2
    y2p = Y[1] - 2 * Y[2] - mp * Y[1] / D1 - mu * Y[1] / D2
    return jnp.array([Y[2], Y[3], y1p, y2p])