Solving partial differential equations (PDEs), equations that model the evolution of physical systems in both space and time, is of ubiquitous interest in applied mathematics and physics. For problems in soft materials, fluctuating hydrodynamics, and more generally, statistical mechanics, solving stochastic partial differential equations (SPDEs) is of particular importance. Although there exists a wide variety of numerical methods to integrate SPDEs in time, many become significantly unstable when taking long time steps (the time step is the length of the time interval over which the state of the system being simulated is updated). Smaller time steps can be used to retain stability and accuracy, but doing so requires more computations to simulate longer stretches of time, which means methods that remain stable for long time steps are desirable for improved efficiency. We utilize the Magnus expansion, which serves as a basis for constructing the solution to a linear differential equation in the form of a matrix exponential, to develop numerical methods for long-time integration of linear SDEs/SPDEs with time-dependent dissipative operators. We introduce efficient methods to compute the random variables in our numerical scheme, ensuring that it adheres to important principles in statistical mechanics including fluctuation-dissipation balance and detailed balance. We compare the performance of our exponential integrator to traditional methods for integrating SPDEs, demonstrating that it remains more stable and accurate for longer time steps. We use our exponential integrator to simulate prototypical dynamical systems including a stochastic pendulum and a fluid subject to shear flow and thermal fluctuations, demonstrating that our methodology can be used to study a broad range of mesoscale phenomena.