@Ritz11 That's a decent answer written from a non-language-specific functional-programming perspective. That is, there are no Maple-specific details other than precisely those copied verbatim from the original Question. So, was that written by ChatGPT or similar AI?

Here it is with some Maple details:

p:= 2*(x-2)^2*(x-1);
F:= factors(p);
seq[fold= (`*`, fun(F[1]))](fun(f[1])^f[2], f= F[2]);
 

@NIMA112 I suspect that the roots are found by trial-and-error. Making implicitplots (as I showed above) helps. Roots are where the red and green cross. You need to change the ranges of A and B. Many choices will show no roots. If you find an intersection, that gives you about 2 significant digits of A and B. You then use these as seed values to fsolve to refine to as many digits as you want (which doesn't always work). The plots are computationally intensive, taking about 2 minutes each.

@sand15 That's strange. I don't know what's happening there. If one were really interested, one could display all the results of the assignments leading up to the creation of the two equations and check where the numeric values in the two executions (2015 and 2023) begin to diverge.

Anyway, this implicitplot clearly shows the 2 roots that I proposed. And by adjusting the ranges appropriately, it'll show the 2 latest roots that @Axel Vogt found. (And it finds another root fairly close to Axel's [A= 37.6324, B= -0.356503].)

plots:-implicitplot(
    Re~(rel), A= -4..1, B= -1..2, 
    gridrefine= 3, crossingrefine= 3, color= [red, green]
);

@sand15 I had moved the OP's Digits:= 50 from where it originally appeared in the worksheet to immediately after the restart: Making that change in your version of the worksheet and re-executing your computation of the residuals using my posted numeric solution yields

eval(rel, {
    A = -2.7553365135418814642586082436429575890825402826031,
    B = -0.70285804987973303586180028708027467941012949957141
});
[-4.426*10^(-43) - 0.*I, 
-4.0506731414438437950164110646723220186973915791808*10^(-40) - 0.*I]

which seems sufficiently close to 0 to me.

@Joe Riel I also found those lines of code, and they clearly show why expand isn't idempotent in this case. The question to me is Shouldn't it be idempotent by design? That is why I only said that perhaps it was a shortcoming or bug. Making it idempotent would be easy.

Yes, I think that the defaults should be changed. Specifically, 

1. It's often difficult to appreciate the 3-dimensional nature of a plot without gridlines, though, unlike your example, I wouldn't necessarily include them on both the inner and outer surfaces unless they were differently colored (which IS fairly easy to do).
2. With the current defaults, the colors of the inner and outer surfaces are too similar to each other.

@rlopez Dr Lopez's suggestion can be given the full sign, primpart, sort treatment to put the results into your requested format, doing them all at once with this single command:

`[]`~(
    L,
    (sign*primpart@sort@lhs@Student:-VectorCalculus:-TangentPlane)~(
        (x - 1)^2 + (y - 3)^2 + (z - 5)^2 = 13^2, `<,>`~(L)
    ) =~ 0
);

I am suprised that this Student command doesn't check whether the specified point is actually on the specified surface. It'll happily return results for points that are not. (This is mentioned on its help page.)

@NIMA112 Why do you suggest this obscure numerical method that'll require much work to get coded into Maple? I'm not saying that it's not a good method for this problem---that I don't know. Do you have any reason to suspect that this method will work better than any of the numerous more-standard methods that are already coded?

@Art Kalb That's a good question, but I don't know the answer. It is perhaps a shortcoming or bug of expand.

@Laurenso Here is the version of the procedure with the on-the-sphere check included:

tangent_plane:= p-> local P:= <p>, R:= <1,3,5> - P;
    [p, `if`(R.R = 13^2, (sign*primpart@sort)(R.(<x,y,z> - P)) = 0, FAIL)]
:

An equivalent to map(expand@expand, N^%H) is 

(expand~@@2)(N^*)

If you want some variables to represent real numbers, then you must use assume or assuming. If you want them all to represent reals, then why not just use the ordinary transpose?

@Laurenso I just realized that my procedure doesn't check whether the given point is actually on the sphere. If you require that check, I can easily add it.

@Laurenso The command to see the procedure is showstat(geom3d::onobjl).

@2cUniverse I think that you've found some relevant info. It'll take me a few days to digest it. Thanks.

Regarding your Reply before the one above: As I said before, that idea is already implemented in the presented code. Once any element x of order k is found, I immediately use it to generate another totient(k) - 1 elements of order k. What I'm missing is a fast way to generate the initial elements x needed to get all elements of order k. Starting with elements of order Carmichael(n) isn't enough.