I try to solve this second order differential equations to get $\phi$:

$$ 4 \Psi \left( \ddot{\phi}_0+ \frac{a^2}{2\phi_0}\dot{\phi}_0 \right) = -\ddot{\phi}+ \left( \frac{a^2 \mathcal{H}}{2\phi_0} \mathcal{H} + \frac{a^2 \dot{\phi}_0}{\phi_0} \right) \dot{\phi}- \left( 3 k^2 - \frac{a^2 \ddot{\phi}_0}{\phi_0} - \frac{\mathcal{H} a^2 \dot{\phi}_0}{\phi_0} - \frac{2 a \dot{a} \dot{\phi}_0}{\phi_0} + \frac{a^2 \dot{\phi}_0}{\phi_0^2}+ \frac{\mathcal{H} a^2}{\phi_0} - \frac{a^4 \dot{\phi_0}}{4\phi^2_0} \right) \phi ~~~~~~~~~~~~ (1) $$

Where $\Psi$ is defined from this ODE:

$$ \left( k^2 + 3 \mathcal{H}^2 +3 \mathcal{H} \frac{a^2\dot{\phi}_0}{2\phi_0} \right) \Psi = - \frac{\ddot{\phi}}{2\phi_0} + \dot{\phi} \left( \frac{\mathcal{H}}{2\phi_0} + \frac{\dot{\phi_0}}{2\phi_0} - \frac{a^2 \mathcal{H}}{2\phi_0} \right) + \phi \left( \frac{\mathcal{H} \dot{\phi_0}}{\phi_0} + \frac{5\mathcal{H}a^2}{2\phi_0}-\frac{a^4 \dot{\phi}_0}{4\phi^2_0} \right) ~~~~~~(2) $$

Both $\phi$ and $\Psi$ are functions of t. $\dot{\phi}=\frac{\partial}{\partial ~t} \phi$, and $\ddot{\phi}=\frac{\partial^2}{\partial^2t} \phi$.

$a(t) = t + \frac{t^2}{2}, H(t) = \frac{\dot{a}}{a},~ \phi_0(t) =- \frac{e^{-Ht}}{H} $ and $k$ is a constant.

So when substituting by (2) into (1) and all the variables, I get a single second-order differential equation in $\phi$

Here is the code that I try:

a[t_]= t + t^2/2; h[t_]= Derivative[1][a][t]/a[t]; f0[t_]= -(Exp[-2*h[t]*t]/h[t]); df0[t_]= D[f0[t], t]; ddf0[t_]= D[f0[t], {t, 2}];

y[t_] := (1/(k^2 + 3*h[t] + 3*h[t]*((a[t]^2*df0[t])/(2*f0[t]))))*(-(Derivative[2][x][ t]/(2*f0[t])) + (Derivative[1][x][t]/(2*f0[t]))*(h[t] + df0[t] - a[t]^2*h[t]) + (x[t]/(2*f0[t]))*(2*h[t]*df0[t] + 5*h[t]*a[t]^2 - (a[t]^4/(2*f0[t]))*df0[t]))

eq := Simplify[-Derivative[2][x][ t] + (Derivative[1][x][t]/f0[t])*((a[t]^2/2)*h[t] + h[t]*f0[t] + a[t]^2*df0[t]) - (x[t]/f0[t])*(3*k^2*f0[t] - a[t]^2*ddf0[t] - h[t]*a[t]^2*df0[t] - 2*a[t]*Derivative[1][a][t]*df0[t] + (a[t]^2*df0[t])/f0[t] + h[t]*a[t]^2 - (a[t]^4*df0[t])/(4*f0[t])) - 4*y[t]*(ddf0[t] + (a[t]^2/(2*f0[t]))*df0[t])]

Where $x(t)=\phi (t)$ , $y(t)=\Psi (t)$ and $f0(t) =\phi_0(t)$. Ofcourse DSolve doesn't work with eq. So I tried:

AsymptoticDSolveValue[{eq == 0, x[0] == 0, Derivative[1][x][0] == 1}, x[t], {t, 0, 10}]

But it also doesn't give a solution. Hint: the initial conditions are arbitrary.

Another Hint: this equation can not be solved numerically because $k$ is an unknown constant.

So I appricticate any help in finding an analytic solution for this equation.