This is a bit messy, but in can be computed using Gröbner basis theory and elimination orders. See Cox, Little, O'Shea "Ideals, Varieties, and Algorithms" for the details on the theory.
We start by writing down some equations:
Suppose we have a point $(a,b)$ and it is tangent to the curve $y = x^2$ at the point $(t,t^2)$. Then
$(t - a)^2 + (t^2 - b)^2 = b^2$ says that the circle and parabola intersect at $t$
The slope of the parabola at $(t,t^2)$ is $2t$ and the slope of the circle at the point $(t, t^2)$ is $-\frac{t-a}{t^2 - b}$ (use implicit differentiation).
Therefore, our two equations are
- $(t - a)^2 + (t^2 - b)^2 - b^2 = 0$ and
- $2t(t^2 - b) + (t - a) = 0$ (it is important to write these equations as "polynomial = 0" for the purpose of feeding them to the computer).
Now, using Macaulay2 (online system here), we have the following session:
i1 : R = QQ[t,a,b]
o1 = R
o1 : PolynomialRing
We are telling Macaulay2 that our variables are t, a and b and that we are working with equations with rational coefficients.
i2 : I = ideal( (t - a)^2 + (t^2 - b)^2 - b^2, 2*t*(t^2 - b) + (t - a) )
o2 : Ideal of R
We tell Macaulay2 that the system of equations will be called "I".
i3 : J = eliminate(t, I)
6 4 2 2 4 4 2 3 4 2 2 2
o3 = ideal(16a - 32a b + 16a b - 40a b - 24a b + a + 12a b - 2a b)
o3 : Ideal of R
We tell Macaulay2 to write the system of equations without using the variable t.
i4 : primaryDecomposition J
2 4 2 2 4 2 3 2 2
o4 = {ideal a , ideal(16a - 32a b + 16b - 40a b - 24b + a + 12b - 2b)}
o4 : List
We tell Macaulay2 to split the ideal into irreducible components.
Macaulay2 has told us that the locus consists of two components. First, is the set $a^2 = 0$ (which gives us all the circles that are tangent at $(0,0)$). The second is the quartic equation
$$ \boxed{16a^4 - 32a^2b^2 + 16b^4 - 40a^2b - 24b^3 + a^2 + 12b^2 - 2b = 0.} $$
Here is what the locus looks like (in black)

See https://www.desmos.com/calculator/p5q7ezcvkg for an interactive version.
Notice that there are two components of the real locus (depending on which side of the circle the tangent intersection happens).
You can also solve for $a$ and $b$ in terms of the intersection parameter $t$. (Same computation in Macaulay2, but eliminate a or b instead of t). This gives
$$ \boxed{a = \frac{1}{2} \left(t\pm\sqrt{4 t^4+t^2}\right), b = \frac{1}{4} \left(4 t^2\pm\sqrt{4 t^2+1}+1\right).} $$
*take $+,+$ and $-,-$ for $t \le 0$ and take $+,-$ and $-,+$ for $t \ge 0$.