Following:
Was the instability with all solvers? Including CVODE and the Fortran methods? Those are the same solvers you can use from Mathematica, so you should get the same results if you choose them.