Set of integrators do not consistently solve an equation in Python

I must solve the following second order differential equation:

$delta phi^{”}_{mathbf{k}}+(3-epsilon)delta phi^{‘}_{mathbf{k}}+left(frac{k^2}{a^2 H^2}+frac{V_{,phiphi}}{H^2}-6epsilon +4epsiloneta -2epsilon^2right)delta phi_{mathbf{k}}=0;$

which is that of scalar field perturbations at linear order in cosmic inflation. $epsilon$, $eta$, $phi$, $V_{,phiphi}$, $H$ and $a$ denote background quantities which can be solved separately (i.e., their corresponding field equations do not contain field perturbation terms). $k$ represents the Fourier mode to which the field perturbation is ascribed and can be viewed as constant.

This equation can also be written as follows:

$frac{1}{a^3 phi^{‘2}H}left[a^3 phi^{‘2} H left(frac{delta phi_{mathbf{k}}}{phi^{‘}}right)^{‘}right]^{‘} =-frac{k^2}{a^2H^2} frac{delta phi_{mathbf{k}}}{phi^{‘}}.$

When I proceed to solve the first equation above, I use $delta phi_{mathbf{k}}$ to solve the second one at the same time. However, it turns out that I do not get to solve the equations consistently (i.e., $a^3 phi^{‘2} H left(frac{delta phi_{mathbf{k}}}{phi^{‘}}right)^{‘}$ obtained from the first equation does not agree with the one got from the last one, by far).

To solve the equations, I am employing various SciPy integrator solvers (RK45, RK23, LSODA, BDF, etc…), and none of them sorts out the discrepancy successfully.

Is there any suggestion about how I could get a consistent (more accurate) solution from those equations? Let me know if further details are needed.