As an example for an ODE course I used the ODE $$ y' = \frac{y}{x} + \frac{1}{\cos(\tfrac{y}{x})} $$ to illustrate domains of existence. Standard substitution $z=y/x$ turns the equation to $$ z' = \frac{1}{x\cos(z)} $$ which can be solved by separation of variables to obtain the solution $$ \phi(x) = x\arcsin(\ln(x/x_0) - \sin(y_0/x_0)) $$ with initial value $\phi(x_0)=y_0$ (and the solution ceases to exists at $x = x_0\exp(1-\sin(y_0/x_0))$).
To check that my solution is correct, I compared the solution with the output of some numerical solver (ode45
from MATLAB). Since I was to lazy to figure out the precise upper bound at which the solution ceases to exist, I just plugged in some upper bound and looked at the result. As expected, something weired happened beyond the point of existence. However, when comparing the numerical solution with the computed one (which is correct, by the way), I noticed something strange: Although the solution computed by the exact formula has a problem (when the argument of the $\arcsin$ becomes larger than one and from that point on MATLAB plot the real part and neglects the imaginary part of the solution), the "wrong exact solution" somehow follows that "wrong numerical solution" in a very crude, strange but not deniable way. Here is the plot for $x_0=1$, $y_0=0$:
(the red is the plot of $\phi$ and the blue is the output of ode45
).
Why is that? What is going on with an ODE solver beyond existence? Why is it following the real part of the non-existing complex solution?