Biggest circle under a polynomial curve

We have

$$ \cases{ y = x(1-x)(2x+1)\\ (x-x_0)^2+(y-y_0)^2=y_0^2 } $$

substituting the first into the second we have

$$ p(x)=(x-x_0)^2+(x(1-x)(2x+1)-y_0)^2-y_0^2 = 0 $$

This polynomial at tangency should have the following structure

$$ p(x) = c(x-x_1)^2(x-x_2)^2((x-a)^2+b^2)\ \ \ \ (\text{Why?}) $$

now comparing coefficients we have

$$ 0=\eta = \left\{ \begin{array}{l} x_0^2-c x_1^2 x_2^2 \left(a^2+b^2\right) \\ 2 c x_1 x_2 \left(x_2 \left(a^2+b^2\right)+x_1 \left(a^2-a x_2+b^2\right)\right)-2 x_0-2 y_0 \\ 2 -c x_1^2 \left(a^2-4 a x_2+b^2+x_2^2\right)-4 c x_2 x_1 \left(a^2-a x_2+b^2\right)-c x_2^2 \left(a^2+b^2\right)-2 y_0 \\ 2 \left(c x_1 \left(a^2-4 a x_2+b^2+x_2^2\right)+c x_2 \left(a^2-a x_2+b^2\right)+c x_1^2 \left(x_2-a\right)+2 y_0+1\right) \\ c \left(a^2+b^2\right)+c \left(4 x_1 \left(x_2-a\right)+x_2 \left(x_2-4 a\right)+x_1^2\right)+3 \\ 2 c \left(x_1+x_2\right)-2 (a c+2) \\ 4-c \\ \end{array} \right. $$

and we finish by solving this nonlinear system for $x_0,y_0,x_1,x_2,a,b,c$ using a minimization procedure such as

$$ \min_{x_0,y_0,x_1,x_2,a,b,c} \eta\cdot\eta\ \ \ \text{s. t.}\ \ \ \{x_0 > 0,y_0 >0, x_1 > 0, x_2 > 0\} $$

Follows a MATHEMATICA script which materializes those computations

Clear[g1, g2]
g1[x_, y_] := y - x (1 - x) (2 x + 1)
g2[x_, y_] := (x - x0)^2 + (y - y0)^2 - y0^2
dif = g2[x, x (1 - x) (2 x + 1)] - c (x - x1)^2 (x - x2)^2 ((x + a)^2 + b^2)
eta = CoefficientList[dif, x]

sol = NMinimize[{eta.eta, x0 > 0, y0 > 0,x1 > 0, x2 > 0}, {x0, y0, x1, x2, a, b, c }]


g2x = g2[x, y] /. sol[[2]]
gr1 = ContourPlot[g2x == 0, {x, 0, 1}, {y, 0, 0.6}, ContourStyle -> Red];
gr2 = ContourPlot[g1[x, y] == 0, {x, 0, 1}, {y, 0, 0.6}];
gr3 = Plot[0, {x, 0, 1}]
Show[gr1, gr2, gr3, AxesOrigin -> {0, 0}, AspectRatio -> 0.6]