I have a system of four coupled nonlinear ODEs. Three of them are of first order, while the fourth equation is of second order. The latter also incorporates a variable coefficient, $a(x)$, which can can take values of $0<a(x)<1$: $$ 0 = a(x) \frac{\partial^2 f(x)}{\partial x^2} + b(x) $$
I'm solving the equation using Matlab. The coefficient is determined from a set of discrete points through interpolation. For most $a(x)$ it works. However, for cases where the coefficient has a value close to 0 or close to 1, the Matlab's solver reports large condition numbers for the Jacobian, or even singular Jacobians, and fails. Manually decreasing the step size in the particular region helps, but not a lot. I don't think it is due to large gradients as the code can handle them as long as the difference in $a$ is not too large.
As this only occurs for extreme values of $a(x)$, I can live with it. Nevertheless, I was wondering what the cause is. Is it that, for $a\approx 1$ the system is of second order and for $a \approx 0$ of first order? I could imagine that this change of order within a short distance causes problems.
I know that stiffness is a very squishy topic, but can anybody give an educated guess or a hint to any literature concerning or mentioning this kind of problem?