A lot of things depend on the properties of $a(t)$, but I will mention some general methods I would use.
For $t \to 0$ the approach is obvious: search for the series solution. If we can represent:
$$x(t)=c_0+c_1 t+c_2 t^2+ \cdots$$
$$a(t)=a_0+a_1 t+a_2 t^2+ \cdots$$
We have:
$$2c_2=c a_0 c_0$$
$$6c_3=c (a_1c_0+ a_0c_1)$$
And so on. With $c_0,c_1$ determined by the initial conditions.
For large times the approach will have to be different. Maybe we can represent $a(t)$ as a series in $1/t$, in which case it will be easy to find the asymptotic solution as a series same as above.
If not, we can try to use various other methods. Like WKB method, mentioned by @LutzL. In fact, you can scroll to an example describing a general equation very similar to your case.
Another approach to try would be using perturbation theory. (We denote $ca(t)=b(t)$ for simplicity.)
Assuming that variable part of $b(t)$ is small, that is, it can be represented as a series in some small parameter $\beta \ll 1$, like this:
$$b(t)=B+\beta b_1(t)+ \beta^2 b_2(t) +\dots$$
$$B=\text{const}$$
Going back to the original equation, we can try to also represent $x(t)$ this way:
$$x(t)=x_0(t)+\beta x_1(t)+ \beta^2 x_2(t) +\dots$$
Then we have:
$$\ddot{x}=b(t) x$$
$$(\ddot{x_0}+\beta \ddot{x_1}+ \beta^2 \ddot{x_2} +\dots)=(B+\beta b_1(t)+ \beta^2 b_2(t) +\dots)(x_0+\beta x_1+ \beta^2 x_2 +\dots)$$
Matching the powers of $\beta$, we have:
$$\ddot{x_0}=B x_0$$
$$\ddot{x_1}=B x_1+b_1(t)x_0(t)$$
$$\ddot{x_2}=B x_2+b_1(t)x_1(t)+b_2(t)x_0(t)$$
And so on. For $x_0$ we obviously have the general solution:
$$x_0(t)=C_1 e^{\sqrt{B} ~t}+C_2 e^{-\sqrt{B} ~t}$$
The other equations all have the same form:
$$\ddot{x_n}=B x_n+f(t)$$
with $f(t)$ a known function. This is a linear inhomogeneous equation and can be solved exactly or approximately by the usual methods. We only need to find a particular solution, because we know the general solution has the same form as $x_0$.
We can use the variable coefficient method, searching for $x_{np}$ as:
$$x_{np}=C(t) e^{\sqrt{B}t}$$
Which brings us to:
$$(\ddot{C}+\sqrt{B} ~\dot{C})e^{\sqrt{B}t}=f(t)$$
$$\dot{C}=D(t)$$
$$\dot{D}+\sqrt{B}~ D=e^{-\sqrt{B}t}f(t)$$
This is a 1st order inhomogeneous equation which is solved in quadratures.
This is a little involved in general, but allows us to find as many orders of approximation as we want.
A different form of the equation might help. Let us change the function:
$$x(t)=e^{y(t)}$$
Then we have:
$$\ddot{y}+\dot{y}^2=b(t)$$
Now introducing another function:
$$z(t)=\dot{y}(t)$$
we have a nonlinear 1st order ODE:
$$\dot{z}+z^2=b(t)$$
This is a Riccati equation. This equation has a number of nice properties, for example if we know any particular solution, we can find the general one.
Riccati equation has also been extensively studied by many researchers for the last 200 years or so, which is a good reason to search online for the methods relevant to your case.