What algorithm/option in DifferentialEquations should I choose to solve efficiently a large and stiff system of ODE (a system of 100 to 10_000 say) when I am *only* interested in the value of the stationary solution (i.e I know that system of ODE converges to a constant vector and I just want to find this constant).

`CVODE_BDF(linear_solver = :GMRES)`

with the solver options `save_everystep=false,save_start=false`

is likely what youâ€™re looking for.

This blog post explains a workflow for doing exactly this: http://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/

If you know that there is a unique stationary solution and you are only interested in it, maybe you try a Newton method to find the equilibrium directly.

True. But that doesnâ€™t capture initial condition dependence if thatâ€™s a thing in this model. FWIW if you use the steady state stuff then this is simplified, and then switching between rootfinding and a dynamical form is trivialized:

http://docs.juliadiffeq.org/latest/solvers/steady_state_solve.html#SteadyStateDiffEq.jl-1

```
# Newton method via nlsolve choice, defaults to NLsolve.jl
sol = solve(prob,SSRootfind())
# Dynamical solve to steady-state, automatically exits when `norm(du)< ...`
sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)
```

DynamicSS is exactly what I was looking for. Thanks a lot!