Thank you Johannes, I've tried your suggestion of changing rtol and atol values but the deSolve solution is always the same. In fact, I have obtained the answer to my question from Thomas Petzoldt and Stephen P. Ellner (thank you!). I copy here the answer from Stephen: "the problem is that this
Cgills=Cgills_ini*exp(-(k01+k12)*t)+(Cu*k01+Crest*k21*wrest/wgills)/(k10+k12)*(1-exp(-(k10+k12)*t)) Crest=Crest_ini*exp(-k21*t)+Cgills*k12*wgills/wrest/k21*(1-exp(-k21*t))
is not the analytic solution to the differential equations. The first line is the solution to the first differential equation, assuming that Crest is constant. The second line is the solution to the second differential equation, assuming that Cgills is constant. But since Cgills and Crest change over time, those are not the solution to the system of differential equations. The answer to trust is the one that you're getting from the ODE solvers in R."
---------------------------------------------------------------------- Message: 2 Date: Thu, 17 Mar 2016 18:26:03 +0100 From: Johannes Ranke<johannes.ranke at jrwb.de> To:r-sig-dynamic-models at r-project.org Cc: Paula Sanchez Marin<paulasanchez at uvigo.es> Subject: Re: [R-sig-dyn-mod] iterative procedure vs deSolve Message-ID: <10725503.32hNF4IJnn at tux> Content-Type: text/plain; charset="iso-8859-1" Dear Paula, did you try to adapt rtol and atol, which are passed from ode() to its default solver lsoda(), to lower values? I had a similar case when I compared an Eigenvalue based solution of a system of ODEs with the solution produced by deSolve, and it turned out the the deSolve solution converged to the Eigenvalue based solution when atol and/or rtol were decreased. Kind regards, Johannes