Over the last decades, factorization theorems became an important method to tackle problems in perturbative QCD, especially within the framework of effective field theories. In Soft-Collinear Effective Theory, these factorization theorems include beam functions accounting for the initial-state collinear interactions. While these functions have been calculated case by case for different observables until now, we are investigating an automated approach for a general class of observables. For this, we study a general phase-space parameterization which factorizes the singularities of the beam function in an universal way. This approach has been implemented in the public code "pySecDec" in order to calculate the next-to-leading order and part of the next-to-next-to leading order beam function.