I don't think there is a solution to those equations, but if there were you could use the follwing method:
import sympy as sp
# Define symbols
theta1, theta2, l1, l2, x, y = sp.symbols("theta1 theta2 l1 l2 x y")
# Define equations, rearranged so expressions equal 0
eq1 = l1 * sp.cos(theta1) + l2 * sp.cos(theta1 + theta2) - x
eq2 = l1 * sp.sin(theta1) + l2 * sp.sin(theta1 + theta2) - y
# Solve for theta1 & theta2
solution = sp.solve([eq1, eq2], [theta1, theta2], dict=True)
print(solution)
I was trying to use the sympy nonlinsolve solver for a similar inverse kinematics problem but noticed this comment in the docs:
Currently nonlinsolve is not properly capable of solving the system of equations having trigonometric functions. solve can be used for such cases (but does not give all solution)