  This paper concerns the analysis of random second order linear differential
equations. Usually, solving these equations consists of computing the first
statistics of the response process, and that task has been an essential goal in
the literature. A more ambitious objective is the computation of the solution
probability density function. We present advances on these two aspects in the
case of general random non-autonomous second order linear differential
equations with analytic data processes. The Fr\"obenius method is employed to
obtain the stochastic solution in the form of a mean square convergent power
series. We demonstrate that the convergence requires the boundedness of the
random input coefficients. Further, the mean square error of the Fr\"obenius
method is proved to decrease exponentially with the number of terms in the
series, although not uniformly in time. Regarding the probability density
function of the solution at a given time, we rely on the law of total
probability to express it in closed-form as an expectation. For the computation
of this expectation, a sequence of approximating density functions is
constructed by reducing the dimensionality of the problem using the truncated
power series of the fundamental set. We prove several theoretical results
regarding the pointwise convergence of the sequence of density functions and
the convergence in total variation. The pointwise convergence turns out to be
exponential under a Lipschitz hypothesis. As the density functions are
expressed in terms of expectations, we propose a symbolic Monte Carlo sampling
algorithm for their estimation. This algorithm is implemented and applied on
several numerical examples designed to illustrate the theoretical findings of
the paper.
