Stiff Stochastic differential equation

Question

I have a stochastic differential equation that has stiffness under some range of parameters. Therefore using numerical simulations for this system can end up with a divergent behavior. Using methods such as BDF and Radau from scipy.optimize.solve_ivp also didn't help with my problem. My guess is that this is because of the stochastic term in my equations. What would you do in such situations? Are there some packages that deal with stiffness and stochasticity at the same time? Or is there a method that can be implemented? Here's the SDE: $$ \begin{aligned} \dot{z_1} = \lambda(t) - z_1 + \sqrt{\lambda}dW ; \dot{z_2} = z_1 - z_2 ; \cdots \dot{z_{n}} = z_{n-1} - z_n \end{aligned} $$ where $$ \lambda(t) = \exp(\beta_1 z_1 + \beta_2 z_2 + \cdots) $$ As you see the $\lambda(t)$ term will make this problem stiff.