Circles inscribed between curves can be specified by a system of equations relative to the coordinates of the center of the circle and the coordinates of the tangent points. Such a system can have 5 or 6 equations and 6 variables, which are mentioned above.
In the case of 5 equations, we can immediately obtain an infinite set of solutions by selecting the ones we need from it. 
(See the attached text for more details.)
The 1st equation is responsible for the belonging of the point of tangency to one of the curves.
The 2nd equation is responsible for the belonging of the point of tangency to another curve.
In the 3rd equation, the points of tangency on the curves belong to the inscribed circle.
In the 4th and 5th equations, the condition is satisfied that the tangents to the curves are perpendicular to the radii of the circle at the points of contact.
The 6th equation serves either to find a specific inscribed circle or to find an infinite set of solutions. It is selected based on the type of curves and their mutual arrangement.

In this example, we search for a subset of the solution set using the Draghilev method by solving the first five equations of the system: we inscribe circles in two "angles" formed by the intersection of the exponent and the ellipse.
The text of this example, its solution in the form of a picture,"big" option and pictures of similar examples.
INSCRIBED_CIRCLES.mw


 

Addition 09/01/24, 
One curve for the first two equations in coordinates x1,x2 and x3,x4
f1:= x1^2 - 2.5*x1*x2 + 3*x2^2 - 1;
f2:= x3^2 - 2.5*x3*x4 + 3*x4^2 - 1;

A simple visual way to show that the parametric equation of a circle is a helix in our three-dimensional space.
Parametric equation of a circle f1 and f2.
The helix is ​​defined by the intersection of two mutually perpendicular cylindrical surfaces f1 and f2.

 

restart; with(plots): 
R := 1.;
f1 := x1-R*cos(x3); 
f2 := x2-R*sin(x3); 
PT := implicitplot3d([f1, f2], x1 = -6 .. 6, x2 = -6 .. 6, x3 = -2 .. 12, numpoints = 10000, style = surface, color = [blue, green], transparency = .5):
IT := intersectplot(f1, f2, x1 = -1 .. 1, x2 = -1 .. 1, x3 = -2 .. 12, thickness = 3, axes = normal, grid = [10, 10, 30]): 
display(PT, view = [-6 .. 6, -6 .. 6, -2 .. 12]);
display(PT, IT, view = [-6 .. 6, -6 .. 6, -2 .. 12]); 
display(IT, view = [-R .. R, -R .. R, -2 .. 12], scaling = constrained)

Does everyone have a good idea of ​​the work of the Draghilev method? For example, in this answer https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions#answer291268
( https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions )
there was a very successful attempt by Rouben Rostamian to calculate the line of intersection of surfaces without applying the Draghilev method.
Let now not 3d, but 8d. And how will the solve command work in this case? Imagine that aij ((i=1..7,j=1..8)) are partial derivatives, and xj (,j=1..8) are derivatives, as in the above example. f8 is responsible for the parametrization condition.
 

restart;
 f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
 f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
 f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
 f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8;
 f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8; 
 f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
 f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;
 f8 := x1^2+x2^2+x3^2+x4^2+x5^2+x6^2+x7^2+x8^2-1; 
allvalues(solve({f1, f2, f3, f4, f5, f6, f7, f8}, {x1, x2, x3, x4, x5, x6, x7, x8}))

And this is how the Draghilev method works in this case.
 

restart; with(LinearAlgebra):
f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8; 
f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8;
f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;

n := 7; 
x := seq(eval(cat('x', i)), i = 1 .. n+1):
 F := [seq(eval(cat('f', i)), i = 1 .. n)]: 
A := Matrix(nops(F), nops(F)+1):
 for j to nops(F) do for i to nops(F)+1 do A[j, i] := op(1, op(i, op(j, F))) 
end do:
           end do: 

# b[i] and b[n+1] are solutions of a linear homogeneous system and at the 
# same time they serve as the right-hand sides of an autonomous ODE.
for i to n do

 b[i] := Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1)) 
                                                                     end do:
 b[n+1] := -Determinant(DeleteColumn(A, n+1)):

Only the original seven linear homogeneous equations with eight variables are needed. We solve them according to Cramer's rule, and in order to have uniqueness when solving the ODE, we use a point on the curve (according to the theory). (This point is sought in any convenient way.)
If we want to get a parameterization, then additionally, directly in dsolve, we can add the following:
 

for i to n do 
b[i] := simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))) 
end do:
b[n+1] := simplify(-Determinant(DeleteColumn(A, n+1))); 
deqs := seq(diff(x[i](s), s) = b[i]/(b[1]^2+b[2]^2+b[3]^2+b[4]^2+b[5]^2+b[6]^2+b[7]^2+b[8]^2)^.5, i = 1 .. n+1):

An example of uniform motion along a generalized coordinate using the Draghilev method. (This post was inspired by school example in one of the forums.)
The equations used in the program are very simple and, I think, do not require any special comments. DM is a procedure that implements the Draghilev method with "partial parameterization".
DM_V.mw

When K = 1, parameterization is carried out by changing the angle of rotation of the wheel. That is, uniform rolling is carried out.

For K = 4, the coordinate corresponding to the position of the slider is parametrized.

When K = 6, the slider moves with acceleration, according to a given equation. Hence, we have carried out the parameterization with respect to “time”.

With the help of such techniques, we can obtain the calculation of the kinematics of both lever mechanisms and various types of manipulators.

It seems to me that Draghilev's method can be applied quite successfully to the solution of Diophantine equations. Here is a simple example where we find two solutions at the intersection line of two ellipsoids:
  x1^2-x1*x2+x2^2+x2*x3+x3^2-961=0;
  (x1-3)^2+10*x2^2+x3^2-900=0;
Solutions: (11, -4, -26) and (10, 1, 29).
 

Based on the text of the program, it is possible to solve various examples with Diophantine equations.
3d_1.mw

Explanations.
f3 is an auxiliary equation for finding the starting point, NPar is a procedure that implements the Draghilev method, the red color of the text is the place where the integer values of the points on the integral curve are filtered.

 Can be compared with the solution of the isolve function
 

 restart:
  f1 := x1^2-x1*x2+x2^2+x2*x3+x3^2-961;
  f2 := (x1-3)^2+10*x2^2+x3^2-900;
  isolve({f1, f2})