Given the cubic Bézier curve of parametric equations:
$$
\begin{aligned}
& b_x(u) := x_0\,(1 - u)^3 + 3\,x_1\,u\,(1 - u)^2 + 3\,x_2\,u^2\,(1 - u) + x_3\,u^3 \\
& b_y(u) := y_0\,(1 - u)^3 + 3\,y_1\,u\,(1 - u)^2 + 3\,y_2\,u^2\,(1 - u) + y_3\,u^3 \\
\end{aligned}
\quad \quad \text{with} \; u \in [0,\,1]
$$
and the circle of parametric equations:
$$
\begin{aligned}
& c_x(v) := x_C + r\,\cos v \\
& c_y(v) := y_C + r\,\sin v \\
\end{aligned}
\quad \quad \text{with} \; v \in [0,\,2\pi)
$$
in order to be tangent it must be:
$$
\begin{cases}
b_x(u) = c_x(v) \\
b_y(u) = c_y(v) \\
b_x'(u) = \lambda\,c_x'(v) \\
b_y'(u) = \lambda\,c_y'(v) \\
\end{cases}
$$
with $0 < u < 1$, $0 < v < 2\pi$, $y_2 \in \mathbb{R}$, $\lambda \in \mathbb{R}$.
In light of this, defined the vector function $\mathbf{f} : \mathbb{R}^4 \to \mathbb{R}^4$ of law:
$$
\mathbf{f}(\mathbf{z}) := \left(b_x(u) - c_x(v), \; b_y(u) - c_y(v), \; b_x'(u) - \lambda\,c_x'(v), \; b_y'(u) - \lambda\,c_y'(v)\right)
$$
with $\mathbf{z} \equiv \left(u,\,v,\,y_2,\,\lambda\right)$, the zeros of this function:
$$
\mathbf{f}(\mathbf{z}) = \mathbf{0}
$$
can be calculated using the Newton-Raphson method, according to which:
$$
\mathbf{z}_{i + 1} = \mathbf{z}_i - \mathbf{J}_{\mathbf{f}}^{-1}(\mathbf{z}_i) \cdot \mathbf{f}(\mathbf{z}_i)
$$
where $\mathbf{z}_0$ is an initial vector sufficiently close to the desired zero; the method stops when:
$$
\left\lVert \mathbf{J}_{\mathbf{f}}^{-1}(\mathbf{z}_i) \cdot \mathbf{f}(\mathbf{z}_i) \right\rVert \le \epsilon\,,
$$
with $\epsilon > 0$ pre-established according to your needs.
Implementing this method (with some tricks) in Wolfram Mathematica 12.2:
{x0, y0} = {0, 10};
{x1, y1} = {0, 1.96};
x2 = 8;
{x3, y3} = {8, 0};
{xC, yC} = {6.713, 5.04};
r = 1.1;
ε = 10^-10;
countmax = 100;
bx[u_] = x0 (1 - u)^3 + 3 x1 u (1 - u)^2 + 3 x2 u^2 (1 - u) + x3 u^3;
by[u_] = y0 (1 - u)^3 + 3 y1 u (1 - u)^2 + 3 y2 u^2 (1 - u) + y3 u^3;
cx[v_] = xC + r Cos[v];
cy[v_] = yC + r Sin[v];
f[{u_, v_, y2_, λ_}] = {bx[u] - cx[v],
by[u] - cy[v],
bx'[u] - λ cx'[v],
by'[u] - λ cy'[v]};
invJf[{u_, v_, y2_, λ_}] = Inverse[Grad[f[{u, v, y2, λ}], {u, v, y2, λ}]];
z = {1/2, 3π/2, (y0 + y1 + y3)/3, 10};
count = 0;
While[Norm[invJf[z].Transpose[{f[z]}]] > ε && count < countmax,
z = z - Transpose[invJf[z].Transpose[{f[z]}]][[1]];
If[z[[1]] < 0 || z[[1]] > 1, z[[1]] = 1/2];
If[z[[2]] < π || z[[2]] > 2π, z[[2]] = 3π/2];
count = count + 1
];
Print[{count, z}]
we get:
{6, {0.689356, 4.18946, 7.66843, 10.7859}}
that is, with only $6$ iterations, we obtain $y_2 = 7.66843$, which is what is desired.
Finally, to verify what has been obtained, by writing:
P0 = {x0, y0};
P1 = {x1, y1};
P2 = {x2, z[[3]]};
P3 = {x3, y3};
P0lab = P0 + {+0.3, +0.2};
P1lab = P1 + {+0.3, -0.2};
P2lab = P2 + {+0.3, +0.2};
P3lab = P3 + {+0.3, -0.2};
plot1 = ParametricPlot[Evaluate[{bx[u], by[u]} /. y2 -> z[[3]]],
{u, 0, 1}, PlotStyle -> Red];
plot2 = ParametricPlot[{cx[v], cy[v]}, {v, 0, 2π}, PlotStyle -> Green];
plot3 = Graphics[{Directive[Blue, Dashed], Line[{P0, P1, P2, P3}]}];
Show[{plot1, plot2, plot3},
AxesLabel -> {"x", "y"}, GridLines -> Automatic, PlotRange -> All,
Epilog -> {Text[Style[Subscript["P", "0"], Blue], P0lab],
Text[Style[Subscript["P", "1"], Blue], P1lab],
Text[Style[Subscript["P", "2"], Blue], P2lab],
Text[Style[Subscript["P", "3"], Blue], P3lab]}
]
we get:

which actually shows what you want. I hope it's clear enough, good luck! ^_^
Addendum 1
A fairly faithful approach to the algorithm implemented in Wolfram Mathematica 12.2 is:
define the parameters $x_0,\,y_0,\,x_1,\,y_1,\,x_2,\,x_3,\,y_3,\,x_C,\,y_C,\,r,\,\epsilon,\,\text{countmax}$;
initialize $\text{count} = 0$, $u = \dots$, $v = \dots$, $y_2 = \dots$, $\lambda = \dots$, then define $\mathbf{z} = \left(u,\,v,\,y_2,\,\lambda\right)^t$;
calculate the components:
$$
\small
\begin{aligned}
& f_1 = x_0 + 3\left(x_1 - x_0\right)u + 3\left(x_2 - 2\,x_1 + x_0\right)u^2 + \left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u^3 - x_C - r\,\cos v \\
& f_2 = y_0 + 3\left(y_1 - y_0\right)u + 3\left(y_2 - 2\,y_1 + y_0\right)u^2 + \left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u^3 - y_C - r\,\sin v \\
& f_3 = 3\left(x_1 - x_0\right) + 6\left(x_2 - 2\,x_1 + x_0\right)u + 3\left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u^2 + \lambda\,r\,\sin v \\
& f_4 = 3\left(y_1 - y_0\right) + 6\left(y_2 - 2\,y_1 + y_0\right)u + 3\left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u^2 - \lambda\,r\,\cos v \\
\end{aligned}
$$
then define $\mathbf{f} = \left(f_1,\,f_2,\,f_3,\,f_4\right)^t$;
calculate the components:
$$
\begin{aligned}
& J_{11} = 3\left(x_1 - x_0\right) + 6\left(x_2 - 2\,x_1 + x_0\right)u + 3\left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u^2 \\
& J_{12} = r\,\sin v \\
& J_{13} = 0 \\
& J_{14} = 0 \\
\\
& J_{21} = 3\left(y_1 - y_0\right) + 6\left(y_2 - 2\,y_1 + y_0\right)u + 3\left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u^2 \\
& J_{22} = -r\,\cos v \\
& J_{23} = 3\,u^2\,(1 - u) \\
& J_{24} = 0 \\
\\
& J_{31} = 6\left(x_2 - 2\,x_1 + x_0\right) + 6\left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u \\
& J_{32} = \lambda\,r\,\cos v \\
& J_{33} = 0 \\
& J_{34} = r\,\sin v \\
\\
& J_{41} = 6\left(y_2 - 2\,y_1 + y_0\right) + 6\left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u \\
& J_{42} = \lambda\,r\,\sin v \\
& J_{43} = 3\,u\,(2 - 3\,u) \\
& J_{44} = -r\,\cos v \\
\end{aligned}
$$
then define $\mathbf{Jf} =
\begin{pmatrix}
J_{11} & J_{12} & J_{13} & J_{14} \\
J_{21} & J_{22} & J_{23} & J_{24} \\
J_{31} & J_{32} & J_{33} & J_{34} \\
J_{41} & J_{42} & J_{43} & J_{44} \\
\end{pmatrix}$;
calculate:
$$
\begin{aligned}
& \mathbf{invJf} = \text{inverse}(\mathbf{Jf})\,; \\
& \mathbf{resf} = \mathbf{invJf}.\mathbf{f} \; \text{[product row by column]}\,; \\
& \text{errf} = \sqrt{\mathbf{resf}_1^2 + \mathbf{resf}_2^2 + \mathbf{resf}_3^2 + \mathbf{resf}_4^2}\,;
\end{aligned}
$$
while $\text{errf} > \epsilon$ and $\text{count} < \text{countmax}$:
$$
\begin{aligned}
& \mathbf{z} = \mathbf{z} - \mathbf{resf}\,; \\
& \left(u,\,v,\,y_2,\,\lambda\right)^t = \mathbf{z}\,; \\
& \text{if} \; u < 0 \; \; \text{or} \; \; u > 1 \; \; \; \text{then} \; u = \frac{1}{2}\,; \\
& \text{if} \; v < \pi \; \; \text{or} \; \; v > 2\pi \; \; \text{then} \; v = \frac{3\pi}{2}\,; \\
& \text{repeat steps 3, 4, 5}\,; \\
& \text{count} = \text{count} + 1\,. \\
\end{aligned}
$$
In particular, to make the calculations faster, it's possible to calculate the residual vector once and for all, reducing everything to a series of elementary operations that are easy to iterate:
(*Input data*)
{x0, y0} = {0, 10};
{x1, y1} = {0, 1.96};
x2 = 8;
{x3, y3} = {8, 0};
{xC, yC} = {6.713, 5.04};
r = 1.1;
ε = 10^-10;
countmax = 100;
(*Residual function definition*)
res[u_, v_, y2_, λ_] := Module[{f1, f2, f3, f4, J11, J12, J21, J22, J23,
J31, J32, J34, J41, J42, J43, J44, num1, num2, num3, num4, den},
f1 = x0 + 3 (x1 - x0) u + 3 (x2 - 2 x1 + x0) u^2 + (x3 - 3 x2 + 3 x1 - x0) u^3 - xC - r Cos[v];
f2 = y0 + 3 (y1 - y0) u + 3 (y2 - 2 y1 + y0) u^2 + (y3 - 3 y2 + 3 y1 - y0) u^3 - yC - r Sin[v];
f3 = 3 (x1 - x0) + 6 (x2 - 2 x1 + x0) u + 3 (x3 - 3 x2 + 3 x1 - x0) u^2 + λ r Sin[v];
f4 = 3 (y1 - y0) + 6 (y2 - 2 y1 + y0) u + 3 (y3 - 3 y2 + 3 y1 - y0) u^2 - λ r Cos[v];
J11 = 3 (x1 - x0) + 6 (x2 - 2 x1 + x0) u + 3 (x3 - 3 x2 + 3 x1 - x0) u^2;
J12 = r Sin[v];
J21 = 3 (y1 - y0) + 6 (y2 - 2 y1 + y0) u + 3 (y3 - 3 y2 + 3 y1 - y0) u^2;
J22 = -r Cos[v];
J23 = 3 u^2 (1 - u);
J31 = 6 (x2 - 2 x1 + x0) + 6 (x3 - 3 x2 + 3 x1 - x0) u;
J32 = λ r Cos[v];
J34 = r Sin[v];
J41 = 6 (y2 - 2 y1 + y0) + 6 (y3 - 3 y2 + 3 y1 - y0) u;
J42 = λ r Sin[v];
J43 = 3 u (2 - 3 u);
J44 = -r Cos[v];
num1 = f1 (J23 (J34 J42 - J32 J44) - J22 J34 J43) +
J12 (f2 J34 J43 + f3 J23 J44 - f4 J23 J34);
num2 = f1 (J23 (J31 J44 - J34 J41) + J21 J34 J43) -
J11 (f2 J34 J43 + f3 J23 J44 - f4 J23 J34);
num3 = (f1 J21 - f2 J11) (J32 J44 - J34 J42) -
(f1 J22 - f2 J12) (J31 J44 - J34 J41) +
(f3 J44 - f4 J34) (J11 J22 - J12 J21);
num4 = (f1 J31 - f3 J11) (J22 J43 - J23 J42) -
(f1 J32 - f3 J12) (J21 J43 - J23 J41) +
(f2 J43 - f4 J23) (J11 J32 - J12 J31);
den = J11 (J23 (J34 J42 - J32 J44) - J22 J34 J43) +
J12 (J23 (J31 J44 - J34 J41) + J21 J34 J43);
Return[{num1, num2, num3, num4} / den]];
(*Attempt solution*)
{u, v, y2, λ} = {1/2, 3π/2, (y0 + y1 + y3)/3, 10};
vres = res[u, v, y2, λ];
errq = vres[[1]]^2 + vres[[2]]^2 +
vres[[3]]^2 + vres[[4]]^2;
count = 0;
(*Update solution*)
While[errq > ε^2 && count < countmax,
{u, v, y2, λ} = {u, v, y2, λ} - vres;
If[u < 0 || u > 1, u = 1/2];
If[v < π || v > 2π, v = 3π/2];
vres = res[u, v, y2, λ];
errq = vres[[1]]^2 + vres[[2]]^2 +
vres[[3]]^2 + vres[[4]]^2;
count = count + 1
];
(*Output data*)
Print[{count, {u, v, y2, λ}}]
Although the definition of these parameters is very tedious, once this work is done, there is the advantage of being able to dominate the algorithm in an optimal way, as there are no more hidden steps.
In particular, I would like to underline the fact that the high efficiency of the Newton-Raphson method is paid for in terms of low robustness, since the convergence of this method is highly influenced by the choice of the first attempt solution vector.
Precisely for this reason, based on the problem under examination, it's necessary to intervene by optimizing it through some more or less intelligent measures. In this case, I intervened inside the while loop with a couple of if conditional checks to make $u$ and $v$ converge to a value within the respective intervals, but this can certainly be done more efficiently, as well as you will almost certainly have to intervene in the testing to adapt it to your needs.
I hope I have explained myself sufficiently, good job! ^_^
Addendum 2
If the specific problem doesn't allow to identify a robust strategy in the choice of the initial vector to trigger the Newton-Raphson method, it's necessary to change direction.
In particular, it's possible to abandon the world of local convergence methods by migrating to that of global convergence methods in which, when it comes to polynomials, the Aberth-Ehrlich method is undoubtedly an excellent compromise between algorithmic simplicity and at the same time efficiency/robustness.
On an algorithmic level, I would proceed like this:
through simple substitutions, the above system of equations can be simplified as follows:
$$
\begin{cases}
a + b\,u + c\,u^2 + d\,u^3 = \cos v \\
e + f\,u + \left(g + 3\,y_2/r\right)u^2 + \left(h - 3\,y_2/r\right)u^3 = \sin v \\
b + 2\,c\,u + 3\,d\,u^2 = -\lambda\,\sin v \\
f + 2\left(g + 3\,y_2/r\right)u + 3\left(h - 3\,y_2/r\right)u^2 = \lambda\,\cos v \\
\end{cases}
$$
with $0 < u < 1$, $0 < v < 2\pi$, $y_2 \in \mathbb{R}$, $\lambda \in \mathbb{R}$;
dividing the first two equations member to member, I obtain the expression of $\sin v/\cos v$, then dividing the last two equations I obtain the expression of $-\sin v/\cos v$ which, when suitably equalized, allow to obtain a quadratic equation in $y_2$, an unknown which can therefore be expressed as a function of $u$;
adding the squares member to member of the first two equations, substituting the expression of $y_2$ and simplifying everything with further squares, it's possible to obtain a polynomial equation $p(u) = 0$, which in the most general case is 12th degree, where the coefficient of the monomial of maximum degree can be made unitary by dividing all the others by itself;
the time has therefore come to bring heavy artillery into the field by applying the above numerical method, which is conceived precisely for monical polynomial equations, in general with complex coefficients, through which it's possible to calculate all twelve complex roots simultaneously by updating them with the following formulation:
$$
u_j' = u_j - \frac{1}{\frac{p'\left(u_j\right)}{p\left(u_j\right)} - \begin{aligned}\sum_{k = 1, \, k \ne j}^{12}\end{aligned} \frac{1}{u_j - u_k}} \;;
$$
all that remains is to examine these roots one by one by selecting only those with a real part between $0$, $1$ and with an imaginary part close to zero, then through simple checks it's possible to determine the set of values $\left(u,\,v,\,y_2,\,\lambda\right)$ that satisfy all four equations placed in the system (due to the way the problem has been formulated, there may be more solutions than there may be, unless we introduce further constraints on the input parameters).
The advantage of this other method is that it doesn't require a particularly reasoned initial vector, it's sufficient to define it once and for all so that the (complex) components are asymmetrical and uniformly distributed (generally assumed on a circle).
On the other hand, it's necessary to refer to the arithmetic of complex numbers (which is the lesser evil, since it's easily algorithmizable), so it's necessary to establish a tolerance with which to discriminate the real roots from the complex roots by comparing it with the imaginary part (which could cause problems, since it isn't so obvious to define it arbitrarily).
Anyway, I tried to write this other algorithm in Visual Basic for Applications (VBA) in Excel and I was very satisfied, it seems very robust, although obviously it can always be improved:
Function sum(x, y)
sum = Array(x(0) + y(0), x(1) + y(1))
End Function
Function diff(x, y)
diff = Array(x(0) - y(0), x(1) - y(1))
End Function
Function prod(x, y)
prod = Array(x(0) * y(0) - x(1) * y(1), x(1) * y(0) + x(0) * y(1))
End Function
Function div(x, y)
div = Array((x(0) * y(0) + x(1) * y(1)) / (y(0) ^ 2 + y(1) ^ 2), (x(1) * y(0) - x(0) * y(1)) / (y(0) ^ 2 + y(1) ^ 2))
End Function
Function pow(x, n)
pow = Array(Cos(n * WorksheetFunction.Atan2(x(0), x(1))) * (x(0) ^ 2 + x(1) ^ 2) ^ (n / 2), Sin(n * WorksheetFunction.Atan2(x(0), x(1))) * (x(0) ^ 2 + x(1) ^ 2) ^ (n / 2))
End Function
Sub solve()
'IMPORT DATA'
x0 = Cells(4, 4)
y0 = Cells(5, 4)
x1 = Cells(7, 4)
y1 = Cells(8, 4)
x2 = Cells(10, 4)
x3 = Cells(13, 4)
y3 = Cells(14, 4)
xC = Cells(16, 4)
yC = Cells(17, 4)
r = Cells(18, 4)
imax = Cells(26, 4)
tol = Cells(27, 4)
segno = Cells(28, 4)
'SCALING'
a = (x0 - xC) / r
b = 3 * (x1 - x0) / r
c = 3 * (x2 - 2 * x1 + x0) / r
d = (x3 - 3 * x2 + 3 * x1 - x0) / r
e = (y0 - yC) / r
f = 3 * (y1 - y0) / r
g = 3 * (-2 * y1 + y0) / r
h = (y3 + 3 * y1 - y0) / r
'POLYNOMIAL COEFFICIENTS'
u0 = 4 * (-1 + a ^ 2) * (-1 + a ^ 2 + e ^ 2)
u1 = -4 * (3 + 3 * a ^ 4 - 3 * a ^ 3 * b + a * b * (3 - 2 * e ^ 2) + e * (-3 * e + f) + a ^ 2 * (-6 + 3 * e ^ 2 - e * f))
u2 = 9 * a ^ 4 + a ^ 3 * (-38 * b + 8 * c) + (-9 + 4 * b ^ 2) * (-1 + e ^ 2) + 8 * a * c * (-1 + e ^ 2) + 14 * e * f - f ^ 2 + a ^ 2 * (-18 + 13 * b ^ 2 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + a * b * (38 + 8 * e * (-3 * e + f))
u3 = 2 * (a ^ 3 * (15 * b + 2 * (-7 * c + d)) + 2 * b * c * (-1 + 2 * e ^ 2) + b ^ 2 * (7 + 2 * e * (-3 * e + f)) + a * (3 * b ^ 3 + 14 * c - 2 * d + b * (-15 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 4 * e * (d * e + c * (-3 * e + f))) + 2 * (f ^ 2 + e * (-3 * f + g + h)) - 2 * a ^ 2 * (b * (11 * b - 4 * c) + f ^ 2 + e * (-3 * f + g + h)))
u4 = -22 * a * b ^ 3 + b ^ 4 + b ^ 2 * (-12 + 37 * a ^ 2 + 10 * a * c + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 2 * (3 * a ^ 3 * (4 * c - 3 * d) + 2 * c ^ 2 * e ^ 2 + a * (c * (-12 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + d * (9 + 4 * e * (-3 * e + f))) - 3 * e * (g + h) + f * (-2 * f + g + h) + a ^ 2 * (2 * c ^ 2 + 2 * f ^ 2 + 3 * e * (g + h) - f * (g + h))) + b * (6 * a ^ 2 * d + 8 * d * e ^ 2 + c * (18 - 62 * a ^ 2 + 8 * e * (-3 * e + f)) - 8 * a * (f ^ 2 + e * (-3 * f + g + h)))
u5 = 2 * (-2 * b ^ 4 + b ^ 3 * c + 9 * a ^ 3 * d + b * (c * (-9 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 4 * d * (1 + e * (-3 * e + f))) + 2 * (c * (d + 2 * d * e ^ 2) + c ^ 2 * (1 + e * (-3 * e + f)) - f * (g + h)) + a ^ 2 * (29 * b * c - 10 * c ^ 2 - 18 * b * d + 2 * f * (g + h)) - 2 * b ^ 2 * (f ^ 2 + e * (-3 * f + g + h)) + a * (10 * b ^ 3 + b ^ 2 * (-22 * c + d) + d * (-9 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 2 * b * (c ^ 2 + 2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)) - 4 * c * (f ^ 2 + e * (-3 * f + g + h))))
u6 = 4 * b ^ 4 - 10 * b ^ 3 * c + 4 * d ^ 2 - 3 * a ^ 2 * d ^ 2 + 4 * d ^ 2 * e ^ 2 + 24 * a * d * e * f - 8 * a * d * f ^ 2 + c ^ 2 * (-6 + 22 * a ^ 2 + 9 * e ^ 2 - 14 * e * f + f ^ 2) - 8 * a * d * e * g - g ^ 2 + a ^ 2 * g ^ 2 - 8 * a * d * e * h - 2 * g * h + 2 * a ^ 2 * g * h - h ^ 2 + a ^ 2 * h ^ 2 + b ^ 2 * (46 * a * c + c ^ 2 - 22 * a * d + 4 * f ^ 2 + 6 * e * (g + h) - 2 * f * (g + h)) - 2 * c * (d * (1 + 9 * a ^ 2 + 12 * e ^ 2 - 4 * e * f) - 6 * a * e * (g + h) + 2 * a * f * (-2 * f + g + h)) + 2 * b * (21 * a ^ 2 * d + d * (-6 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + a * (-13 * c ^ 2 - 2 * c * d + 4 * f * (g + h)) - 4 * c * (f ^ 2 + e * (-3 * f + g + h)))
u7 = 2 * (b ^ 3 * (6 * c - 2 * d) - 3 * c * d + 15 * a ^ 2 * c * d - 3 * d ^ 2 + 9 * c * d * e ^ 2 - 6 * d ^ 2 * e ^ 2 + 6 * c ^ 2 * e * f - 14 * c * d * e * f + 2 * d ^ 2 * e * f - 2 * c ^ 2 * f ^ 2 + c * d * f ^ 2 - 2 * c ^ 2 * e * g - 2 * c ^ 2 * e * h + b ^ 2 * (16 * a * d - c * (4 * c + d) + 2 * f * (g + h)) - 2 * a * (c ^ 3 + c ^ 2 * d - 3 * d * e * (g + h) - 2 * c * f * (g + h) + d * f * (-2 * f + g + h)) + b * (2 * c * (2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)) + a * (17 * c ^ 2 - 8 * c * d - 3 * d ^ 2 + (g + h) ^ 2) - 4 * d * (f ^ 2 + e * (-3 * f + g + h))))
u8 = 8 * b ^ 3 * d + 9 * a ^ 2 * d ^ 2 + 9 * d ^ 2 * e ^ 2 + 24 * c * d * e * f - 14 * d ^ 2 * e * f + 4 * c ^ 2 * f ^ 2 - 8 * c * d * f ^ 2 + d ^ 2 * f ^ 2 + 6 * c ^ 2 * e * g - 8 * c * d * e * g - 2 * c ^ 2 * f * g - 2 * c * (4 * d * e + c * (-3 * e + f)) * h + b ^ 2 * (13 * c ^ 2 - 2 * c * d - 2 * d ^ 2 + (g + h) ^ 2) + 2 * b * (-c ^ 3 - c ^ 2 * d + d * (3 * a * d + 4 * f ^ 2 + 6 * e * (g + h) - 2 * f * (g + h)) + c * (22 * a * d + 4 * f * (g + h))) + 2 * a * (4 * c ^ 3 + c ^ 2 * d + 4 * d * f * (g + h) + c * (-3 * d ^ 2 + (g + h) ^ 2))
u9 = 2 * (2 * b ^ 2 * d * (4 * c + d) + 2 * c ^ 2 * f * (g + h) + 2 * c * d * (2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)) + a * d * ((7 * c - d) * (c + d) + (g + h) ^ 2) + b * (3 * c ^ 3 + 2 * c ^ 2 * d + 6 * a * d ^ 2 - c * d ^ 2 + 4 * d * f * (g + h) + c * (g + h) ^ 2) - 2 * d ^ 2 * (f ^ 2 + e * (-3 * f + g + h)))
u10 = c ^ 4 + 2 * c ^ 3 * d + 2 * c * d * (3 * a * d + 5 * b * d + 4 * f * (g + h)) + c ^ 2 * (10 * b * d + d ^ 2 + (g + h) ^ 2) + 2 * d * (2 * b ^ 2 * d + b * (g + h) ^ 2 + d * (3 * a * d + 2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)))
u11 = 2 * d * (c ^ 3 + 2 * c ^ 2 * d + 2 * d * (b * d + f * (g + h)) + c * (2 * b * d + d ^ 2 + (g + h) ^ 2))
u12 = d ^ 2 * ((c + d) ^ 2 + (g + h) ^ 2)
'VECTORING'
ReDim coeff(1 To 1, 1 To 13)
coeff(1, 1) = Array(u0, 0)
coeff(1, 2) = Array(u1, 0)
coeff(1, 3) = Array(u2, 0)
coeff(1, 4) = Array(u3, 0)
coeff(1, 5) = Array(u4, 0)
coeff(1, 6) = Array(u5, 0)
coeff(1, 7) = Array(u6, 0)
coeff(1, 8) = Array(u7, 0)
coeff(1, 9) = Array(u8, 0)
coeff(1, 10) = Array(u9, 0)
coeff(1, 11) = Array(u10, 0)
coeff(1, 12) = Array(u11, 0)
coeff(1, 13) = Array(u12, 0)
'COEFFICIENTS ANALYSIS'
n = 13
flag = 0
While flag = 0
If Abs(coeff(1, n)(0)) < 10 ^ -15 Then
n = n - 1
Else
flag = 1
End If
Wend
'NORMALIZE COEFFICIENTS'
For i = 1 To n - 1
coeff(1, i) = div(coeff(1, i), coeff(1, n))
Next i
n = n - 1
'ATTEMPT u-SOLUTION'
ReDim sol(1 To 1, 1 To n)
For i = 1 To n
sol(1, i) = Array(Cos(2 * WorksheetFunction.Pi() * i / n), Sin(2 * WorksheetFunction.Pi() * i / n))
Next i
'UPDATE u-SOLUTION'
For i = 1 To imax
'ABERTH-EHRLICH METHOD'
For j = 1 To n
'GENERATE p'
p = Array(0, 0)
For k = 1 To n
p = sum(p, prod(coeff(1, k), pow(sol(1, j), k - 1)))
Next k
p = sum(p, pow(sol(1, j), n))
'GENERATE dp'
dp = Array(0, 0)
For k = 2 To n
dp = sum(dp, prod(Array(k - 1, 0), prod(coeff(1, k), pow(sol(1, j), k - 2))))
Next k
dp = sum(dp, prod(Array(n, 0), pow(sol(1, j), n - 1)))
'GENERATE q'
q = Array(0, 0)
For k = 1 To n
If k <> j Then
q = sum(q, div(Array(1, 0), diff(sol(1, j), sol(1, k))))
End If
Next k
sol(1, j) = diff(sol(1, j), div(p, diff(dp, prod(p, q))))
Next j
Next i
'SEARCH y2-SOLUTION'
i = 1
flag = 0
While flag = 0 And i <= n
If sol(1, i)(0) > 0 And sol(1, i)(0) < 1 And Abs(sol(1, i)(1)) < tol Then
u = sol(1, i)(0)
a0 = 9 * u ^ 3 * (2 - 3 * u) * (1 - u)
b0 = 3 * u * ((2 - 3 * u) * (e + f * u + g * u ^ 2 + h * u ^ 3) + u * (1 - u) * (f + 2 * g * u + 3 * h * u ^ 2))
c0 = (a + b * u + c * u ^ 2 + d * u ^ 3) * (b + 2 * c * u + 3 * d * u ^ 2) + (e + f * u + g * u ^ 2 + h * u ^ 3) * (f + 2 * g * u + 3 * h * u ^ 2)
cosine = a + b * u + c * u ^ 2 + d * u ^ 3
If b0 ^ 2 - 4 * a0 * c0 > 0 And Abs(2 * a0) > 10 ^ -15 And cosine > -1 And cosine < 1 Then
y2 = r * (-b0 - (b0 ^ 2 - 4 * a0 * c0) ^ (1 / 2)) / (2 * a0)
sine = e + f * u + (g + 3 * y2 / r) * u ^ 2 + (h - 3 * y2 / r) * u ^ 3
v = (1 - segno) * WorksheetFunction.Pi() + segno * WorksheetFunction.Acos(cosine)
If Abs(sine - Sin(v)) < tol Then
flag = 1
lambda = -(b + 2 * c * u + 3 * d * u ^ 2) / sine
End If
If flag = 0 Then
y2 = r * (-b0 + (b0 ^ 2 - 4 * a0 * c0) ^ (1 / 2)) / (2 * a0)
sine = e + f * u + (g + 3 * y2 / r) * u ^ 2 + (h - 3 * y2 / r) * u ^ 3
If Abs(sine - Sin(v)) < tol Then
flag = 1
lambda = -(b + 2 * c * u + 3 * d * u ^ 2) / sine
End If
End If
End If
End If
i = i + 1
Wend
'PRINT SOLUTION'
If flag = 1 Then
Cells(4, 21) = u
Cells(5, 21) = v
Cells(11, 4) = y2
Cells(6, 21) = y2
Cells(7, 21) = lambda
Else
Cells(4, 21) = ""
Cells(5, 21) = ""
Cells(11, 4) = ""
Cells(6, 21) = ""
Cells(7, 21) = ""
Output = MsgBox("Solution not found!", vbCritical, "WARNING!")
End If
'OPTIMIZE ASPECT RATIO'
xmin = Application.WorksheetFunction.Min(Union(Range("D4:D4"), Range("D7:D7"), Range("D10:D10"), Range("D13:D13"), Range("F4:F104"), Range("I4:I104")))
xmax = Application.WorksheetFunction.Max(Union(Range("D4:D4"), Range("D7:D7"), Range("D10:D10"), Range("D13:D13"), Range("F4:F104"), Range("I4:I104")))
ymin = Application.WorksheetFunction.Min(Union(Range("D5:D5"), Range("D8:D8"), Range("D11:D11"), Range("D14:D14"), Range("G4:G104"), Range("J4:J104")))
ymax = Application.WorksheetFunction.Max(Union(Range("D5:D5"), Range("D8:D8"), Range("D11:D11"), Range("D14:D14"), Range("G4:G104"), Range("J4:J104")))
minabs = Application.WorksheetFunction.Min(xmin, ymin) * 110 / 100
maxabs = Application.WorksheetFunction.Max(xmax, ymax) * 110 / 100
ActiveSheet.ChartObjects("Grafico 1").Activate
ratio = 5 / 6
If xmax - xmin > ratio * (ymax - ymin) Then
ActiveChart.Axes(xlCategory).MinimumScale = minabs
ActiveChart.Axes(xlCategory).MaximumScale = maxabs
ActiveChart.Axes(xlValue).MinimumScale = minabs / ratio
ActiveChart.Axes(xlValue).MaximumScale = maxabs / ratio
Else
ActiveChart.Axes(xlCategory).MinimumScale = minabs * ratio
ActiveChart.Axes(xlCategory).MaximumScale = maxabs * ratio
ActiveChart.Axes(xlValue).MinimumScale = minabs
ActiveChart.Axes(xlValue).MaximumScale = maxabs
End If
ActiveSheet.Range("A1").Select
End Sub

I hope I haven't written nonsense, in case let me know in the comments. Good luck. ^_^