I'm not sure about the problem you want to solve. The system of ODEs has a unique solution, given the initial conditions. So, do you want to find the initial conditions such that all constraints are satisfied? In that case, you can probably make some headway analytically. The first constraint implies that $\dot{q}_1 = \dot{q}_2$ and $\ddot{q}_1 = \ddot{q}_2$ so $f_1(q_1,q_2) = f_2(q_1,q_2)$; that should simplify your problem.
Or, perhaps you already know that the constraints are satisfied, because they are a consequence of the differential equation (and initial conditions), and you want to make sure that the numerical solution also satisfies the constraint. That is studied in "Geometric Integration", and a good book is Hairer, Lubich & Wanner, Geometric Numerical Integration, Springer. The projection method that Sasha mentioned is one possibility, but in fact many methods (like Runge-Kutta) follow linear constraints automatically.
Or, perhaps your description is not complete, and you have for instance three unknown functions, described by two differential equations and one algebraic constraint. Such problems are called differential-algebraic equations (DAEs). One book on their numerical solution that I can recommend is: Ascher & Petzold, Computer methods for ordinary differential equations and differential-algebraic equations, SIAM, 1998.