A two-point boundary value problem whose highest-order derivative is a Riemann–Liouville–Caputo derivative of order \(\alpha \in (1,2)\) is considered. A similar problem was considered in Gracia et al. (BIT 60:411–439, 2020) but under a simplifying assumption that excluded singular solutions. In the present paper, this assumption is not imposed; furthermore, the finite difference method of the BIT paper, which was proved to attain 1st-order convergence under a sign restriction on the convective term, is replaced by a piecewise polynomial collocation method which can give any desired integer order of convergence on a suitably graded mesh. An error analysis of the collocation method is given which removes the above sign restriction and numerical results are presented to support our theoretical conclusions. The tools devised for this analysis include new comparison principles for Caputo initial-value problems and weakly singular Volterra integral equations that are of independent interest. Numerical experiments demonstrate the sharpness of our theoretical results.
