Needs to be replaced with:
After that, note that the equation for dy(1) evaluates to NaN at f = 0.5 and to -Inf at f = 0 (I only explored between 0 and 1). In the original question it was stated that the the tspan input to ode45 is [0.5 1]. So that will cause a problem because the initial value of dy is NaN at f = 0.5 and that kills the entire solution. I ran from [0.51 1] and got a non-NaN result. In the follow up comment the user enters the tspan limits. If those limits cause the solver to hit either f = 0 or f = 0.5 exactly there will be a problem. So I guess the real question is if it's expected that the equations for dy have these singulariites and how they should be handled if the solver has to integrate through them.