@one man Yes, it's about precision.
It is easy to exhibit a circle with r = 0.7698144869481772563348519  (Digits=25) having a unique tangency point. But visually is already obvious.

sol := [0.75, [r = 0.75, x0 = 1.08, x1 = 0.699729170630851, x2 = 1.27446720511737, x3 = 1.89574141571275, y0 = 1.16, y1 = 0.505162888466810, y2 = 1.89574141571275, y3 = 1.27446720511737]]:
p1:=plots:-implicitplot(e, 0..2, 0..2): p2:=plots:-implicitplot(h, 0..2, 0..2):
p3:=plottools:-circle(eval([x0,y0],sol[2]), eval(r,sol[2]), color=red):
p4:=plots:-pointplot(eval([[x0,y0],[x1,y1],[x2,y2],[x3,y3]],sol[2]), symbolsize=8, color=blue):
p5:=plots:-textplot(eval({seq}([x||i,y||i,P__||i], i=0..3), sol[2]), align={below, left}, color=blue):
plots:-display(p1,p2,p3,p4,p5);

(I modified sol manually. I have not computed, but visually it seems I am right.)

@one man  If you impose a unique point of tangency on each curve, then there is no circle with maximal radius. 
In fact, for any r' < r = 0.7698... (and r' > 0.75, say) we can decrease the radius of the circle to r'  and move it a bit to be tangent just near P₁ and near P₂ .

@nm I don't think so. You simply apply it outside your context.

Let's say you have a function F defined on sets, but F({g,{h}})  is not necessarily equal to F({g,h}).
Then, if you apply Carl's FlattenSet to

r:={a,{b,c},d,{e,f,F({g,{h}})}};

you will come up here too with "This code has serious bug in it".

This is a problem in number theory, not quite a Diophantine equation; an IMO level one, so not very simple.
It can be transformed into a Diophantine equation:

isolve( x^2+y^2 = k*(1+x*y));

          {k = _Z1^2, x = _Z1, y = 0}

Actually this is given by Maple 2018 (on an old Laptop of mine). This is not correct without x>=0, y>=0,  for example x=2, y=-1, k=-5.
In Maple 2024 the result is NULL (also wrong, of course).

Note that isolve is not even able to solve directly the equation for k=4,  although it is of a Pell type. But we can help here:

isolve(z^2-3*y^2 - 4 = 0);  # x = 2*y + z

  ==>  z= ..., y=...

So, unfortunately isolve is not useful for our problem. The only thing we can do is to check numerically the result:

for i to 1000 do for j to i-1 do
  if irem(  ( i^2+j^2) , 1+i*j)=0 then lprint([i,j],(i^2+j^2)/(1+i*j)) fi
od od:

[8, 2], 4
[27, 3], 9
[30, 8], 4
[64, 4], 16
[112, 30], 4
[125, 5], 25
[216, 6], 36
[240, 27], 9
[343, 7], 49
[418, 112], 4
[512, 8], 64
[729, 9], 81
[1000, 10], 100

@mmcdara You have chosen the HermiteH polynomials in the Galerkin method. They are orthogonal in the interval (-oo, oo) with the weight exp(-x^2) you did not use.
So, why HermiteH? Probably ChebyshevT with the corresponding weight are better.

Most  coulditbe(f(z)=0) assuming real;  for not very simple f(z) are anyway (and inevitable) FAIL.

@acer  If we must start (for some reason) with Digits = 15, then method=_d01ajc  will not work, and the nonNAG methods are really slow. 

@sija  Your problem is to solve equations of the form  Int( f(t), t=-4 .. x) = a, f(t) and a being given.

The function F(x) = Int( f(t), t=-4 .. x)  is the solution of the ode  F'(t) = f(t), F(-4) = 0. So, F can be obtained by dsolve (numerically very efficient).

Finally, fsolve(F - a, -10..10)  finds the wanted approximate solution (in the interval -10..10)).
 

@Exiu 

convert(f, radical);

@Christian Wolinski 

coulditbe(exp(z)=0);        #? true
coulditbe(abs(z)+2+z=0);    #? true
solve(abs(z)+2+z=0);        #ok NULL

f can be represented this way iff 0 is not a pole for f (i.e. 0 is not a root of denom(f)).

@Andiguys  I have already mentionrd that Maximize does not accept double inequalities a <= b <=c; use a<=b, b <= c.
So, correct C2:=....

@Andiguys  You do not have a constraint plane. Using
solve(C1 union C3 union C3);
it results
dom := tau1=0..lambda, lambda= 0.0224489798..1.022448980;

Then:
plot3d(TRC(tau1, lambda), dom, orientation=[165, 75,0]);

@Andiguys  Use:

plot3d(TRC(tau1, lambda), tau1 = 0 ..1, lambda = 0 .. 1, orientation=[165, 75,0]);