Getting complex solutions with MATLAB's ode45 I'm trying to solve the ODE below using Matlab's ode45:
$$\frac{dy}{dx}=\sqrt{A\cdot \sin(y) + B\cdot \cos(y)}$$
where $0\le x \le 1$.
The Matlab code is:
A = -2.45;
B = 2.50;
tspan  = 0:0.01:1;
odefun = @(t,y) sqrt(A*sin(y) + B*cos(y));
[t,y] = ode45(odefun, [0 1], 0);

The constants $A$ and $B$ are chosen such that inside the square root it is always non-negative. The solution seems fine for most values of $x$, but for the last few points close to $x=1$ I get complex $y$ with very small imaginary parts. I then took the real or the absolute value of $y$, but then I would get negative value inside the square root.
Should I try other ODE solvers? Or is there anything else I could do to eliminate the complex solutions?
 A: If you square the equation and take the derivative you get
$$
2y'y''=(A\cos(y)-B\sin(y))y'
$$
so that on segments where $y$ is not constant you get the regular ODE
$$
y''=\frac12(A\cos(y)-B\sin(y)),~~y(0)=0,~~y'(0)=\sqrt{B}
$$
which is easily solved numerically

For the solution of the original equation, one would have to detect the maxima/minima and continue from there as a constant solution. 

These connection points of the composite solution only have a continuous first derivative, which invalidates the assumptions of higher-order solvers, and more so of step size controllers.
More specific to this problem, the evaluations of the stages take place close by but not necessarily on the solution curve. In this case this can lead to a negative value for the term under the root, even if this term for the exact solution is non-negative.
Another cause for unexpected behavior could be that the method errors accumulate so that the numerical solution does not stay below the critical value. Here in the solution above with standard tolerances this happens in the interval 
$(0.953,0.975)$ where $A\sin(y) + B\cos(y)$ takes values down to $-0.0003$. So try to set sharper tolerances.
A: This first-order ODE is such that $dy/dx \geqslant 0$. Therefore, a solution $y$ must be an increasing function of $x$. Here, we start with the initial condition $y(0)=0$, where the derivative satisfies $dy/dx \approx 1.6 > 0$. The derivative $dy/dx$ becomes complex when its square becomes negative. For this to happen, an equilibrium value (i.e. where the derivative vanishes) must have been reached. Here, the smallest positive equilibrium value is
\begin{equation}
y^* = {-\arctan(B/A)} \approx 0.796 \, .
\end{equation}
This can be observed when using Matlab: the value $y^*$ is reached numerically just before $x=1$. The choice $A=2.45$ and $B=2.5$ does not avoid $dy/dx$ to become complex.
Note: the Picard-Lindelöf theorem applies here. If $x^*$ is such that $y(x^*)=y^*$, we know that the solution $y$ over $[0, x^*[$ is unique. Thus, it can be approximated numerically with classical methods such as ode45.
A: In reviewing my previous answer it became apparent that the solution there is not the most simple one.
The most straight-forward solution is to extend the domain of the differential equation without changing it on the given domain, so that reaching and slightly crossing the boundary of the given domain does not lead to non-real values for the derivatives and then the values of the numerical solution. This in the main part is achievable here by replacing the square root with an extension to the whole real line. For the solution to stabilize at the boundary one can take the odd extension
$$
r(u)={\rm sign}(u)\sqrt{|u|}.
$$
Numerical solvers find large or infinite values for derivatives difficult to deal with, leading to very small step sizes. To mollify this one can replace some small segment around $u=0$ with a line like in
$$
r(u)=\frac{u}{\max(\epsilon,\sqrt{|u|})}
$$
or similar constructs. To get large step sizes in the asymptotic segment of the solution one will need implicit solvers. Using the implicit "Radau" solver gives the following picture where the internal steps are noted with crosses:

