syms c1 c2 k1 k2 l1 l2 x1 x2
eqns = [k1*sin(c1)+k2*sin(c2-c1) == l1*cos(x1) - l2*cos(x1+x2)
k2*cos(c2-c1)-k1*cos(c1) == l1*sin(x1) - l2*sin(x1+x2)]
sol = solve(eqns, [x1 x2])
This may take some time when it is first done -- about 90 seconds.
The result sol.x1 is two atan() expressions that differ in sign at a single location about 1600 characters into the 4000-ish character expression.
The result sol.x2 is two atan() expression that differ in leading sign (that is, one is the negative of the other)
In the form of solve() above you will get a warning
Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'
In practice, if you ask for return conditions, it will think for a fair while... and then tell you that it was not able to come up with any solutions at all.
That should be taken as a hint to verify the results numerically a couple of times before relying on the solution; it might be the case that these are false roots.
