Indeed it works if I set Digits:=20 for stiff=true or rosenbrock so I can now use range and store the values. One question. I will extend this system to many more coupled ODEs. Say about 10000. However, I am interested in plotting only the first few, or combinations of the first few f[n,l]s, such as here f[0,0] and f[1,0]. In other words I wont need stored all the solutions for f[n,l] with arbitrary n and l in the range xini to xmax. With the option range the solution is stored, at least that is my understanding (all the f[n,l]s are computed and stored, so I might reach a problem with memory). So, how to I keep most of the dsolve solutions as procedures that are estimated once we need a their values for a certain x and only a few of them stored in memory for fast accessing when needed, such as the f[0,0](x), f[0,1](x) in the example above.