For some reason, I would need to know what are the timelike geodesic congruences which cover the entire exterior region of the Schwarschild spacetime. In fact the only thing I really need is the geodesic vector fields defining these congruences.

I have tried something, but being an absolute non-expert in this field I'm not quite sure about it. Here is my attempt:

The derivative with respect to an affine parameter of a generic timelike geodesic is

$$\left\{\matrix{\dot t&=&(1-{r_s\over r})^{-1}E\cr
\dot r&=&\pm \left(E^2-(1-{r_s\over r})(1+{L^2+Q\over r^2})\right)^{1/2}\cr
\dot\theta&=&\pm\left(Q-L^2\mathrm{cotan}^{2}\theta\right)^{1/2}r^{-2}\cr
\dot\varphi&=&{L\over r^2\sin^2\theta}
}\right.$$

where $E$, $Q$, $L$ are constants. Seeing the RHS as a vector field, we obtain a geodesic congruence of each value of the triple of constants. However, we must have $L=0$ for the field to be defined everywhere in the outside region. Hence  the solution to our problem is given by the familly of vector fields above with $L=0$ and $E,Q$ arbitrary.

Is this correct or am I missing something ?

Aside: The texts I have found concerning the Schwarzschild geodesics always focus on the equatorial ones for which $\theta$ is constantly equal to $\pi/2$. I found the formula above in the Wikipedia entry for the Kerr metric, and I set $a$ to $0$. I would  appreciate any reference for this formula in the context of the Schwarzschild metric.