We consider the problem of constructing solutions to the fractional Yamabe problem that are singular at a given smooth sub-manifold, for which we establish the classical gluing method of Mazzeo and Pacard ([63]) for the scalar curvature in the fractional setting. This proof is based on the analysis of the model linearized operator, which amounts to the study of a fractional order ODE, and thus our main contribution here is the development of new methods coming from conformal geometry and scattering theory for the study of non-local ODEs. Note, however, that no traditional phase-plane analysis is available here. Instead, first, we provide a rigorous construction of radial fast-decaying solutions by a blow-up argument and a bifurcation method. Second, we use conformal geometry to rewrite this non-local ODE, giving a hint of what a non-local phase-plane analysis should be. Third, for the linear theory, we use complex analysis and some non-Euclidean harmonic analysis to examine a fractional Schrödinger equation with a Hardy type critical potential. We construct its Green's function, deduce Fredholm properties, and analyze its asymptotics at the singular points in the spirit of Frobenius method. Surprisingly enough, a fractional linear ODE may still have a two-dimensional kernel as in the second order case.