# Solving first order ODE by Runge-Kutta method

I want to solve the first order ODE by the Runge-Kutta method of order 4 $$y'(x)=1+(E-V(x)-1)\sin^2(y(x)),\qquad{}V(x):=\begin{cases} -20\left(1+36x^2\right)(x+1)^2(x-1)^2&\hbox{if }|x|<1,\\ 0&\hbox{else.} \end{cases}$$ with initial condition $y(0)=\pi/2$ where $E$ (eigenvalue of the Hamiltonian with potential $V$) is constant.

So far I have tried to program the Runge-Kutta method of order 4 as

rkstep[f_, y_, h_] := Module[{r1, r2, r3, r4},
r1 = h f[y];
r2 = h f[y + r1/2];
r3 = h f[y + r2/2];
r4 = h f[y + r3];
y + (r1 + 2 r2 + 2 r3 + r4)/6]

and produce a list of function values at each step for a fixed step size $x/nn$

rk[f_, x_, y0_, nn_] :=
NestList[rkstep[f, #, x/nn] &, y0, nn]

Now, I need to specify the ODE by means of a function which I can apply the Runga-Kutte method to - and that's where I get stuck.

So, how do I define an apporpriate function $f$ for the ODE which I can plug in to the Runge-Kutta method defined above?

Edit: I have defined (and omitted the energy)

f[{y_,xx_}] := {1 + (-V /. x -> xx - 1)*(Sin[y])^2, 1}

Then, the following code

oddsol[n_] := {#[[2]], #[[1]]} & /@ ruku[f, {Pi/2, 0.}, 1, n]
ListLinePlot[oddsol[100], PlotRange -> All]

seems to give the desired approximation. However, I need to check for which energies the final condition $$y(1)=\arctan\left(-\frac{1}{\sqrt{-E}}\right)+n\pi$$ is met. But I have no clue, how to implement a parameter energy into the Runge-Kutta method.

-
NDSolve[] is a built-in function. If you absolutely must use fourth-order RK, see this or this. –  Ｊ. Ｍ. Jun 14 at 17:24
That is where I have adapted the code from. However, I cannot figure out how to define the ODE for fourth-order RK. In particular I have problems with the x-dependency of my potential. –  gofvonx Jun 14 at 17:42
The usual trick is to let the eigenvalue be a separate function, using the knowledge that it's constant. Since you know what the derivative of a constant is, it should be easy to write the second differential equation you need. BTW, you seem to have not yet specified the particular potential function you're interested in. –  Ｊ. Ｍ. Jun 14 at 17:45
I still don't understand why you don't want to use NDSolve[], to be honest... –  Ｊ. Ｍ. Jun 16 at 14:40
Because I absolutely must use fourth-order RK. ;-) –  gofvonx Jun 17 at 17:37