This paper deals with random fractional differential equations of the form, CD0+αX(t)+AX˙(t)+BX(t)=0, t>0, with initial conditions, X(0)=C0 and X˙(0)=C1, where CD0+αX(t) stands for the Caputo fractional derivative of X(t). We consider the case that the fractional differentiation order is 1<α<2. For the sake of generality, we further assume that C0, C1, A and B are random variables satisfying certain mild hypotheses. Then, we first construct a solution stochastic process, via a generalized power series, which is mean square convergent for all t>0. Secondly, we provide explicit approximations of the expectation and variance functions of the solution. To complete the random analysis and from this latter key information, we take advantage of the Principle of Maximum Entropy to calculate approximations of the first probability density function of the solution. All the theoretical findings are illustrated via numerical experiments.