Using Matlab, I am trying to solve numerically the following nonlinear system of ODEs:
$$\begin{aligned} \dot B &= -\alpha B -\nu BV \\ \dot X &= A-\mu_1 X -c E(B)VX \\ \dot Y &= -\mu_2 Y +c E(B)VX \\ \dot V &= kY - (\gamma B + E_0)V - \delta V \\ \end{aligned}$$
with
$$E(B)=(\gamma B +E_0)e^{-\beta \gamma B}$$
The system tries to capture the dynamics of how the body of an infant reacts to the presence of the dengue virus.
$V$ is the number of virus.
$B$ the number of antibody.
$X$ the number of healthy cells.
$Y$ the number of infected cells.
$E$ is a function enhancing the effect of antibodies, depending on the number of antibodies.
Implementing this model and trying different solver, I noticed that solvers ode15s and ode23s are performing way better than ode45. Hence, I deduced that my problem was stiff. But with certain set of parameters I had the following warning from Matlab:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.334107e-017.
I worried for a while but I found the malicious ill-conditionned matrix. It's actually that solvers ode23s and ode15s use the Jacobian of the system. I have then computed myself this Jacobian and found the following:
$$DJ(B,X,Y,V)= \begin{pmatrix} -\alpha - \nu V & 0 & 0 & - \nu B \\ -c V X E'(B) & -c V E(B)-\mu_1 & 0 & -c X E(B) \\ c V X E'(B) & c V E(B) & -\mu_2 & c E(B) X\\ -\gamma V & 0 & k & -(\gamma B+E_0+\delta) \\ \end{pmatrix} $$
with
$$ E'(B)=e^{-\beta \gamma B}(\gamma-\beta \gamma^2 B-E_0 \beta \gamma)$$
at this point you have to know a bit more about the dynamics of the system. There is two stable points , one where the virus dies , the antibody disappears , the cells stays healthy
$$P_0=(0, \frac{A}{\mu_1},0,0)$$
And another one, depending of the pararameters in a more complex way :
$$ P_1 = \left(0,\mu_2 \frac{E_0+\delta}{k c E_0},\bar{Y}, \frac{k}{E_0+\delta} \bar{Y }\right) \text{ with } \bar{Y}= \frac{A}{\mu_2} - \frac{\mu_1(E_0+\delta)}{k c E_0} $$
which means, basically, that some cells become infected and stay infected, the virus population explodes before stabilising, and the antibodies disappear.
This second point is the problem because lots of initial conditions seems to be attracted by it and it means that $V$ becomes huge (around $10^{12}$) when the total number of cells doesn't exceed $10^{7}$ and the population of antibodies stays below $10^{4}$.
When $V$ becomes huge, I check the Jacobian again, knowing that $\gamma \approx 1$. We have this one value in the left down corner becoming really huge compared to the other values, also the parameter $c$ is around $10^{-10}$, so really it's mostly about this left corner term.
From there I'm not sure what to do. Is there a way to make a time-varying, well-conditioned Jacobian matrix? Or do you know about any other numerical method that could fit my needs?
It doesn't have to be implemented in Matlab. For now, the results I get are mainly garbage.