But the method I use is more direct, and I think mathematically more robust. By direct I mean it calculates the result directly and does not depend on intermediate results, such as calculating the line or one of x and y first. As such there's much less chance the order of calculation will affect the result.
The approach relies on the fact that either point of intersection forms a triangle with the two circle centres (OTP in the diagram, from Wikipedia) And the three sides of this triangle are known, as they are just the two radii OT and TP with the distance between the centres OP.
The area of a triangle can be calculated from its sides using Heron's formula. Because the area is ½(base × height), and the base OP the known, the height TQ can be calculated. This splits the triangle into two right angled triangles, e.g. OTQ. So using Pythagoras's theorem the length OQ can be calculated.
As a further check the length QP can be calculated the same way, and it can be confirmed their sum is OP. But given the lengths OQ and TQ the coordinates of the point of intersection can be calculated, as the OQ × the normalised vector between the centres + QT × the Perp() of this vector. For the other intersection point change the '+' to a '−'.