You're missing a term in your variation of parameters formula. Given an nth order linear ODE
$\displaystyle a_n(x)y^{(n)} + a_{n-1}(x)y^{(n-1)} + \cdots + a_1(x)y' + a_0(x)y = f(x)$
with independent solutions $\displaystyle y_1, \cdots , y_n$ the variation of parameters formula is
$\displaystyle \displaystyle y_p = \sum_{k=1}^n y_k \int \dfrac{W_k}{W} \dfrac{f(x)}{a_n(x)}dx$.
With what you have you set $\displaystyle a_n = 1$ where it should be $\displaystyle a_n = x^3$. With this assignment things work out giving you the desired $\displaystyle y_p$.