I’m using Mathematica home edition software to numerically solve a specific inflation equation in cosmology. The ODE equation is forth- order, non-linear, stiff. I was using the stiffness switching method and with Working Accuracy of 60, 100, 200, but I’m not getting a stable result. Plugging in the result back into the equations, I’m getting a large deviation, especially around the area that I wish to extract the info. It is a computation issue, and I wonder if someone can recommend different method, or a different software to try. I'm sure there are plenty of researchers that use to deal with similar cosmology inflation equations all the time. thx Ezra
This is not an answer to your question, but more of an observation: More often than not, an ODE is "stiff" or a linear system is "nearly singular" because of a mistake either in deriving the equation to be solved, or in implementing it. Trying to find a way to solve what you have is then just a way to paper over the problem. If you were to find a way to stably solve for what you're interested in, you'd get a number that is, nevertheless, meaningless.
As a consequence, I would recommend to go back to the derivation of where your model comes from and to really understand why it is stiff, and if that is something that you expect for physical reasons. Only when you are convinced that the equation should indeed be stiff does it make sense to think about what method you should use to solve it.
Empirically, in 90% of the cases I've seen, the problem is then in fact no longer stiff/illconditioned (as this was an artifact/bug) and in the other 10% thinking long enough about the origin of the problem suggests how one should approach solving the problem.