Simulating mixed-state evolution in open quantum systems is crucial for various chemical physics, quantum optics, and computer science applications. These simulations typically follow the Lindblad master equation dynamics. An alternative approach known as quantum state diffusion unraveling is based on the trajectories of pure states generated by random wave functions, which evolve according to a nonlinear Itô–Schrödinger equation (ISE). This study introduces weak first-order and second-order solvers for the ISE based on directly applying the Itô–Taylor expansion with exact derivatives in the interaction picture. We tested the method on free and driven Morse oscillators coupled to a thermal environment and found that both orders allowed practical estimation with a few dozen iterations. The variance was relatively small compared to the linear unraveling and did not grow with time. The second-order solver delivers a much higher accuracy and stability with bigger time steps than the first-order scheme, with a small additional workload. However, the second-order algorithm has quadratic complexity with the number of Lindblad operators as opposed to the linear complexity of the first-order algorithm.