Add separate solver for inhomogeneous term that during numerical integration is equal to zero

[clinssen]

If the dynamics of the system is defined such that there is an inhomogeneous (constant-term) component in an ODE, but if numerically, if the inhomogeneous part is zero (this could happen in the case that one or more parameters get set to zero), then because of the way the update expressions are formulated, you could get a division by zero and this results in a NaNs in the output.

To fix this, ODE-toolbox should automatically check whether setting the inhomogenous part to zero results in a division by zero, and if this is the case, return a separate solver. This is akin to checking for singularities in the propagators, but it occurs for the update expressions instead.

The new condition format could look like this (for instance, for a parameter $p$ in the inhomogeneous part):

{
    "solver": "analytical",
    "state_variables": [
        "z"
    ],
    "initial_values": {
        "z": "1"
    },
    "condition": {
        "p==0": {
            "propagators": {
                "__P__z__z": "[...]"
            },
            "update_expressions": {
                "z": "[...]"
            }
        },
        "otherwise": {
            "propagators": {
                "__P__z__z": "[...]"
            },
            "update_expressions": {
                "z": "[...]"
            }
        }
    }
}