import sympy as sy
q0,q1,x,y = sy.symbols('q_0,q_1,x,y')
s = sy.Matrix([sy.Eq(x-sy.cos(q0)-sy.cos(q0+q1),0),
            sy.Eq(y-sy.sin(q0)-sy.sin(q0+q1),0)]);
sol = sy.solve(s.subs({x:2,y:0}),q0,q1)