I am not an experienced Mathematica user, but decided to use it this time because of the in-buit special functions that it has. 

My problem is a solution of the following ODE

    y''[x]=-f(y[x]),  y'[0]=y'[L]=c,

with $L>0$, $c>0$, and periodic rhs $f(y+2\pi)=f(y)=-\Im\{\text{Li}_{2}(-e^{iy})\}$. I tried to solve the equation with the in-built `NDSolve` function and encountered the following problem. 

Using the first integral of the above equation we find that

    0.5*(y'[x])^2+F(y[x])=E,

where $E$ is a constant, and $F(y)=F(y+2\pi)$ such that $F'(y)=f(y)$. Using the boundary conditions we infer that

    F(y[0])=F(y[L]),

implying that `y[L]=y[0]+2*Pi*n`. While solving the equation numerically with `NDSolve`, this property is not respected. In particular 

    s[L_, c_] := 
      NDSolve[{\[Phi]''[y] == -Im[PolyLog[2, -E^(I \[Phi][y])]], \[Phi]'[
      0] == c, \[Phi]'[L] == c}, \[Phi], {y, 0, L}];
    a = s[5, 0.2];
    Evaluate[\[Phi][5] /. a] - Evaluate[\[Phi][0] /. a]

gives me an output $5.02913$, which is not nearly a multiple of $2\pi$. My guess is that something is wrong with the quadrature. Can you help me please? Maybe I have to use a special method, like a shooting method, or something? If so then what is the deal, are there any caveats?

Best,
Kiryl