In this paper we study the randomized heat equation with homogeneous boundary conditions. The diffusion coeffcient is assumed to be a random variable and the initial condition is treated as a stochastic process. The solution of this randomized partial differential equation problem is a stochastic process, which is given by a random series obtained via the classical method of separation of variables. Any stochastic process is determined by its finite-dimensional joint distributions. In this paper, the goal is to obtain approximations to the probability density function of the solution (the first finite-dimensional distributions) under mild conditions. Since the solution is expressed as a random series, we perform approximations of its probability density function. We use two approaches: broadly speaking, first, dealing with the random Fourier coefficients of the random series, and second, taking advantage of the Karhunen-Loeve expansion of the initial condition stochastic process. Finally, several numerical examples illustrating the potentiality of our findings with regard to both approaches are presented.