CompleteSquare(x^2+y^2-2*x-y-2 = 10, x);

restart:
EQ1:=-1958143.922*k*wr+2468.8339*k^3*wr-0.9481118254e16*k^2-114000.8376*k^4:
EQ2 :=-1186578.220*R*k^2*wr-312683.0293*k^5-288960.9621*k^3*R:
 allvalues(solve([EQ1, EQ2], [k, wr]));

Looks like he's alive and well. Those who have lost hope of communicating with this person can find him, for example, here and here.

Look, suddenly this will help.
https://www.maplesoft.com/applications/view.aspx?SID=149514

allvalues(solve({eq1, eq2}, {A, B}));

For calculating the kinematics of the manipulator, Maple of almost any version is quite suitable directly. If a system of nonlinear equations is used to describe the model of a manipulator, then it is very easy mathematically to fix any desired degree of freedom.
Perhaps you can be useful the messages, which can be found at this link:
https://www.mapleprimes.com/posts/210003-Manipulator-With-Variable-Length-Of

 

restart;
f := d+(c*x^3+b*x^2)*(x-1)+(b*x^2+d)^3*(a+x);
expand(f);
sort(collect(expand(f), x), x, ascending);
sort(collect(expand(f), x), x, descending);

There is a universal approach in plots[implicitplot3d]. It can be used also for other equations.
(In the text, any combination of multiplication of two equations.)
composite_surfaces_MP.mw

 implicitplot(x^2*y^2-2*x*y^3+y^4-y^3+x^2 = 1, x = 100 .. 1000, y = 0 .. 1000, numpoints = 100000);

The first that came to mind

restart: with(plots): with(plottools):
c := fieldplot3d([r*sin(t)*cos(p), r*sin(t)*sin(p), r*cos(t)], r = 1 .. 2, t = 0 .. 2*Pi, p = 0 .. 2*Pi, coords = spherical, grid = [25, 25, 25]):
s := sphere([0, 0, 0], 1, color = white, style = surface):
display(c, view = [-1 .. 1, -1 .. 1, -1 .. 1]);
display(c, s, view = [-1 .. 1, -1 .. 1, -1 .. 1])

>? trig
combine(A*B,trig);

For example:
simplify(sin(x)*cos(y)+cos(x)*sin(y),trig);
and
combine(sin(x)*cos(y)+cos(x)*sin(y),trig);
 

It seems very simple:

eps=0.001;
lp:=dsolve(....,type = numeric, events = [[x(t) = 0-eps, [halt]]], range = ta .. te);

(and search the forum, for example: "halt")

For example, nops(solve(x^2+x+1,[x]));

nops(solve([f1(x1,...,xn),..., fn(x1,...,xn)],[x1,...,xn]));

To find the roots of the system of equations NxN, the method works as follows. The initial values of the system of equations NxN we define ourselves. An additional variable with initial value 1 is introduced,  as in the homotopy method, and after that, solving a known system of differential equations, we follow the change of sign of this additional variable. The place where the new variable changes sign be the solution of the original system of equations NxN. 
sistema9_1.mw
As we see, several solutions were found for one initial approximation.

As for the normalization of variables in the differential equation, there is no need for this way of applying  the Diaghilev  method.

And another example with two equations
https://www.mapleprimes.com/posts/200585-Draghilevs-Method-Fx0-Animation