We develop a fourth-order simulation algorithm for solving the stochastic Langevin equation. The method consists of identifying solvable operators in the Fokker-Planck equation, factorizing the evolution operator for small time steps to fourth order, and implementing the factorization process numerically. A key contribution of this paper is to show how certain double commutators in the factorization process can be simulated in practice. The method is general, applicable to the multivariable case, and systematic, with known procedures for doing fourth-order factorizations. The fourth-order convergence of the resulting algorithm allowed very large time steps to be used. In simulating the Brownian dynamics of 121 Yukawa particles in two dimensions, the converged result of a first-order algorithm can be obtained by using time steps 50 times as large. To further demonstrate the versatility of our method, we derive two new classes of fourth-order algorithms for solving the simpler Kramers equation without requiring the derivative of the force. The convergence of many fourth-order algorithms for solving this equation are compared.
