ode45 does not simply compute the function at "every" point and do something like trapazoidal rule. Instead, it starts from a particular set of conditions, and does careful probing around with different times and different initial conditions, hoping the calculated error tolerences come out low enough that it considers that it has been successful in projecting over a range of values. If the error it calculates comes out small enough, it accepts the projected time and boundary conditions and works forward from there. If, though, the error it calculates comes out too large, then it declares the step failed, and goes back to where it was and takes a smaller step. As steps are successful, it tends to grow the step size, and as steps fail, it shrinks the step size. In this way, it can find rapid changes in slope and be more careful in those cases; however, this algorithm can also fail to notice a sufficiently narrow peak if the step size happens to be large enough at the time and the value on the other side of the peak is similar enough to what it is expecting.
The locations it samples at to decide whether to accept or reject are according to mixing coefficients carefully calculated from theory -- there is no randomness or arbitrary choices involved other than the initial step size. It uses the theoretically best coefficients, according to certain assumptions about continuity. There is no room for Mathworks to "tune" the coefficients for choosing where to project to; these are the best locations.
However, since it is doing projections and cannot know ahead of time whether any one projection is reasonable (since the ode function is a "black box" to it, not something it can examine the mathematics of), on projections that are going to end up failing, it can end up asking about the curve values at combinations of boundary conditions that are in locations steeply off of the real curve -- locations that you wouldn't think to ask about if you were taking smaller steps, but which are reasonable to ask about when you realize that the steps are designed to process the curve as rapidly as is compatible with being able to keep reasonable error bounds and this is just probing hoping this direction might be a good one. Projection is successful most of the time, but sometimes you can be close enough to boundary conditions that the function value has to suddenly be large.
Your function happens to have an exp(-y(1)./es), and if y is large enough and negative, then the -y(1)/es can be large and positive and the exp() can go to infinity, leading to an infinite uv result. Because exp() is non-linear, it does not necessarily take much curve projection to go from "hiking along smoothly" to "fell off of the cliff!" if you take too ambitious of a step.
When you generate the inf, then it gets cycled into forming the next set of input coefficients, typically resulting in inf in several input coefficients, and it then being common to end up subtracting or dividing infinities, and both of those lead to nan. And once you have a nan, everything goes out the window, "nan poisoning".
Now, in such cases where you overflow, if somehow you stopped short of infinity, but with a sufficiently large finite value, then the calculation of error would proceed, and MATLAB would see the error was too big and would fail the step, and reduce the stepsize and try again, and the step that would have produced the infinity would be thrown away as irrelevant. The infinity is not a problem because it is infinite: it is a problem because it typically leads to multiple infinities that lead to calculations between infinities that give nan, and nan is a problem. If the function produced an infinity at just the right point or if ode45() somehow noticed it in time, then the step would just be rejected and the algorithm would be back up and try again as designed.
So, what you can do is detect the infinity yourself and return a large enough value to "poison" the calculation to be reasonably sure that the error generated would be too large and that the code will back up. In the case of your function, you can do that by adding the lines:
if any(~isfinite(uv))
uv = [1e10;1e10] .* sign(uv) ;
end
This keeps the result finite but large enough to be over error tolerances (usually.) This line of code is not fool-proof though; in other equations it might have to be adjusted for different error values.
So.. detect the infinity, send out a large value to create a situation where the error will be large enough that ode45 will back up, but without having sent out an actual infinity that can lead to infinity mixes that lead to nan poisoning.
Best Answer