Maybe this could be of help? http://lh3lh3.users.sourceforge.net/solveode.shtml

Generally, speaking LSODA (ODEPACK) is widely used by a number of languages, including R and Python and it can solve both stiff and non-stiff, deciding itself which way to go; the reason that it makes it so popular to my view. Most of the languages have their  own wrappers for the FORTRAN library (much faster implementation), which means that depending on the language you are using, you must find the relevant libraries to load. For example, for Python you may want to read this http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.odeint.html or for R this http://www.inside-r.org/packages/cran/deSolve/docs/lsoda. Other implementations/libraries may be available. 

Hope it helps.

Hi there,

I used Sundials with C/C++ and it was ok.

https://computation.llnl.gov/casc/sundials/main.html

Regards,

I would recommend to consult the book of Hairer and Wanner: Solving Ordinary Differential Equations II: Stiff and Differential Algebraic Problemss.

https://books.google.de/books?id=m7c8nNLPwaIC&pg=PR4&lpg=PR4&dq=hairer+wanner+norsett&source=bl&ots=-93xQHJDAh&sig=o3MaHydvxjDvL_CFwNT8nLZ2Jrs&hl=de&sa=X&ei=X-urVJGEN-TOyQO8uYH4DA&ved=0CFAQ6AEwBzgK#v=onepage&q=hairer%20wanner%20norsett&f=false

Ernst Hairer will surely be able to tell you what might be the best solver for your problem in a prescribed language.

Kind regards,

Hans Josef

Not a bad idea is to use the MATLAB, but you have to select a proper function for stiff equations. MATLAB is a very nice tool for ODE.

An even better idea  is to use Mathematica, extremely stiff equations can easily be solved using it, without to much additional work.

Hello,

you may also consider the RADAU5 method which a (fully implicit) Runge-Kutta method see http://www.unige.ch/~hairer/software.html

Best

ode15s in matlab is widely used  since it is adaptive and well suited for stiff problems

You can use Gear's method for stiff systems of ODEs. Check the IMSL package for an implementation. Or you can use the fundamental matrix method that is adapted to stiff systems. You can refer to the paper

O. Asfar and A. Hussein, "Numerical solution of linear two-point boundary problems via the fundamental matrix method", Int. J. Num. Meth. Engr., vol. 28, No. 5, pp. 1205-1216, May 1989,.

Although this paper is concerned with two-point problems, you can adapt it to initial-value problems as we did in the following paper:

M. Bataineh and O. Asfar, "Application of Multiple Scales Analysis and the Fundamental Matrix Method to the Design of Apodized Rugate Filters:  Initial-Value and Two-Point Boundary Problem Formulations",  IEEE Journal of Lightwave Technology , vol. 18, No. 12, pp. 2217-2223, Dec. 2000.

Besides MATLAB ode15s, you may also use one of MATLAB  ode23s, ode23tb.

These solvers are alright when one is dealing with pure ODEs stiff models. The situation gets complicated when PDEs are reduced to ODEs using method of lines approach in the case of Neuman Boundary conditions.

This restricts applicability of these softwares.

I recommend  exponential integrators for large stiff ODE systems. These integrators are based on computing exponential matrices using krylov subspace methods.  They  are quite popular in numerical analysis community.    Please read  the abstract of the following   standard reference  

http://epubs.siam.org/doi/abs/10.1137/S1064827595295337.

But   dozens of papers are published after this  paper.   if it is suitable for your problem , you may find implementation of the method in some  enviroments.

May be you should reconsider your discretization method. if  all well-known stiff  solvers fail, it is possible for your matrix  involved in computations to be  quite ill conditioned.  if it is so,  none of the solvers works.

thank you for your clear  explanation .

see the paper  below by trefethen and kassam on exponential integrators with matlab codes enabling you do exercises  about your problem. 

fourth order time stepping for stiff pdes, siam journal on scientific computing

may be it is not so relevant but i wonder

from what type of pde do you obtain your ode system ? and  which space integration method    do you use to get  your ode system.    for   example if you have a hyperbolic conservation law, the combination of weno in space and tvd runge kutta in time performs quite well.  a literature review of concerning  your pde may provide deeper insight.

I recommend the book Hall G., Watt J.M. "Modern Numerical Methods for Ordinary Differential Equations".

Article G. Hall and J. M. Watt (editors), Modern Numerical Methods f...

FortranCalculus, http://fortrancalculus.info/apps/fc-compiler.html , has many solvers that may help. It's freeware! Try it and see if it does the job you need to solve.

Two main integration strategies are today adopted for solving stiff ODE systems: the variable-coefficient Backward Differentiation Formula (BDF) methods and the extrapolation methods.

BDF methods were introduced in computational chemistry since the 1970’s and they are still today the most popular method. In BDF methods, a first-order approximation with a very small time step is adopted at the beginning of the integration. Then, as the integration proceeds, a weighted average of information from previous internal steps is used to increase the order of the time-integration scheme ad to increase the time step. Both integration order and time steps can increase significantly during the integration. In particular, if the integration is allowed to continue until steady-state conditions, BDF methods are extremely efficient. In BDF methods, the computational effort is concentrated at the beginning of the integration. The most popular open source BDF solvers on the web are listed in the following:

• CVODE (https://computation.llnl.gov/projects/sundials/cvode): it is a solver for stiff and non-stiff problems. The methods used in CVODE are variable-order, variable-step multistep methods. For non-stiff problems, CVODE includes the Adams-Moulton formulas, with the order varying between 1 and 12. For stiff problems, CVODE includes the BDFs in so-called fixed-leading coefficient form, with order varying between 1 and 5. For either choice of formula, the resulting nonlinear system is solved (approximately) at each integration step;
• DASPK (http://www.cs.ucsb.edu/~cse/software.html): it is a solver for systems of differential-algebraic equations. It includes options for both direct and iterative (Krylov) methods for the solution of the linear systems arising at each (implicit) time step;
• DVODE (https://computation.llnl.gov/casc/odepack/): it is a general purpose solver very similar to DLSODE (see below). However, it uses variable-coefficient methods (fixed-leading coefficient form) instead of the fixed-step-interpolate methods in DLSODE. This and other features make it often more efficient than DLSODE;
• DLSODE (https://computation.llnl.gov/casc/odepack/): it solves stiff and non-stiff systems. Adams methods (predictor-corrector) are adopted in the non-stiff case, and the Backward Differentiation Formulas (BDF) methods (the Gear methods) in the stiff case. The linear systems that arise are solved by direct methods (LU factorization/solution). DLSODE superseded the older GEAR and GEARB packages, with some algorithmic improvements;
• DLSODA (https://computation.llnl.gov/casc/odepack/): it is variant version of the DLSODE package. It switches automatically between stiff and non-stiff methods. This means that the user does not have to determine whether the problem is stiff or not, and the solver will automatically choose the most appropriate method. It always starts with the non-stiff method;
• OdeSMOKE++ (http://cpc.cs.qub.ac.uk/summaries/AEVY_v1_0.html): it is a solver for very stiff problems, based on the BDF in the so-called fixed-leading coefficient form, with order varying between 1 and 5. For either choice of formula, the resulting nonlinear system is solved (approximately) at each integration step. The OdeSMOKE++ solver is part of the OpenSMOKE++ framework is based on the same numerical approach adopted by the BzzMath library;
• RADAU5 (http://archimede.dm.uniba.it/~testset/solvers/radau5.php): it uses an implicit Runge-Kutta method of order 5 (three stages) with step size control and continuous output.

In extrapolation methods, each time step is successively portioned in sub-intervals, a low-order method is as adopted over each sub-interval, and a recursive procedure is used to project to higher-order accuracy over each. Partitioning continues until the user-specified level of accuracy is reached. Similar to BDF methods, the next time step is estimated during the integration on the basis of theoretical and heuristic arguments. In extrapolation methods the computational effort is not concentrated at the beginning of the time interval and in order to advance from time level n to n+1, only information from level n is needed. One of the most popular extrapolation method is SEULEX (http://www.unige.ch/~hairer/prog/stiff/seulex.f).