We present a new algorithm for computing the quasi-stationary distribution of subcritical Galton-Watson branching processes. This algorithm is based on a particular discretization of a well-known functional equation that characterizes the quasi-stationary distribution of these processes. We provide a theoretical analysis of the approximate low-rank structure that stems from this discretization, and we extend the procedure to multitype branching processes. We use numerical examples to demonstrate that our algorithm is both more accurate and more efficient than other approaches.