Computing the stochastic entropy production associated with the evolution of a stochastic dynamical system is a well-established problem. In a small number of cases such as the Ornstein-Uhlenbeck process, of which we give a complete exposition, the distribution of entropy production can be obtained analytically, but in general it is much harder. A recent development in solving the Fokker-Planck equation, in which the solution is written as a product of positive functions, enables the distribution to be obtained approximately, with the assistance of simple numerical techniques. Using examples in one and higher dimension, we demonstrate how such a framework is very convenient for the computation of stochastic entropy production in diffusion processes.