I tried your system and it seems to be a problem with the definition of the system itself. First, I tried applying an L-stable integrator which should be stable on any problem, especially at low tolerances
sol = solve(prob,Rosenbrock23(),reltol=1e-12,abstol=1e-12)
With that failing, I checked the suspicion that it could just be our algorithms by using Sundials’ CVODE stiff integrator:
sol = solve(prob,CVODE_BDF(),reltol=1e-12,abstol=1e-12)
Which also fails. Hairer’s methods also fail. So everybody’s implementation fails at around the same time point, indicating an issue with the problem’s formulation. To check that it wasn’t numerical issues, I used very high precision calculations and still had the equation go unstable at the same point:
u0 = big.([10^-5,1])
tspan = (10^4,10^-4)
p=[0.2]
prob = ODEProblem(FlowEqs,u0,tspan,p)
sol = solve(prob,Rodas5(),reltol=1e-18,abstol=1e-18)
For reference, CVODE_BDF is the CVODE method from Sundials, the C++ library and Hairer’s rodas, radau, etc. are the Fortran methods employed in most other libraries. So this tests not just the Julia-based methods but also the standards everywhere else, with the same results. This is strong empirical information that the derivative function that has been given defines not just a stiff ODE but an unstable ODE. Without knowing the problem itself I am not sure how the ODE should be changed, but hopefully that solves any doubt you may have.