# computation of an integral for 2nd order non degenerate perturbation theory

I am given that the potential of a diatomic molecule is equal to $$V(rho)=-2V left ( frac{1}{rho ^2}-frac{1}{2 rho ^2} right )$$ With $$rho=r/a$$ is a dimesionless coordinate, and $$r$$ is the separation distance between the two atoms. I found the first order corrections without issue, but I am stuck on finding the second order one. I know that
$$E_n^2=sum_{mneq n}frac{|langlepsi_m^0|H’|psi_n^0rangle|^2}{E_n^0-E_m^0}$$
For the given problem, I also found that the wavefunctions are $$psi=e^{-x^2/2} H_{n}(x)$$ (hermite polynomials). My problem is in computing the integral of the inner product, i.e, computing
$$langlepsi_m^0|H’|psi_n^0rangle=int_{-infty}^{infty} x^3 e^{-x^2} H_{n}(x)H_{m}(x)dx$$
I could apply integration by parts a bazillion times to obtain the answer, but it is far too tedious. From reading Griffiths, I know that there is a way to do this much simpler with dirac notation and ladder operators, in conjunction with the usual ladder operator identities, but I am unsure on how to go about this.