The command  Get('xxx1')   returns a string, but you need an exprseq. Use parse(Get('xxx1')).
Or, simply add the line
xxx := parse(xxx);  o := parse(o);  co := parse(co);

These are mathematical conventions. For degree(0) alternative conventions are -1  and "undefined".
- infinity is preferable because e.g. the equality  degree(f*g) = degree(f) + degree(g)  remains valid.
These conventions are related (and compatible) with other two standard ones:
max({}) = - infinity  and  min({}) = infinity.

Try the old-fashioned solution:

f:=x->x^2:
printf("a=?");
a:=parse(readline());
'f(a)'=f(a);

Using a translation, we may suppose that the center of the sphere is (0,0,0).
Even so, the number of the solutions is large and not very easy to obtain (the brute force is out of the question).
E.g.  for the radius of the sphere in the interval  2.10 .. 2.11  here are 32 solutions (actually almost 6 times greater because the permutations are not included; probably other symmetries should have been excluded too):

{{[-20, -3, -3], [-20, 3, 3], [20, -3, 3], [20, 3, -3]}, {[-20, -3, 3], [-20, 3, -3], [20, -3, -3], [20, 3, 3]}, {[-19, -3, -3], [-19, 3, 3], [19, -3, 3], [19, 3, -3]}, {[-19, -3, 3], [-19, 3, -3], [19, -3, -3], [19, 3, 3]}, {[-18, -3, -3], [-18, 3, 3], [18, -3, 3], [18, 3, -3]}, {[-18, -3, 3], [-18, 3, -3], [18, -3, -3], [18, 3, 3]}, {[-17, -3, -3], [-17, 3, 3], [17, -3, 3], [17, 3, -3]}, {[-17, -3, 3], [-17, 3, -3], [17, -3, -3], [17, 3, 3]}, {[-16, -3, -3], [-16, 3, 3], [16, -3, 3], [16, 3, -3]}, {[-16, -3, 3], [-16, 3, -3], [16, -3, -3], [16, 3, 3]}, {[-15, -7, -1], [-15, 7, 1], [15, -5, -5], [15, 5, 5]}, {[-15, -7, -1], [-5, -13, -9], [5, 15, 5], [15, 5, 5]}, {[-15, -7, 1], [-15, 7, -1], [15, -5, 5], [15, 5, -5]}, {[-15, -7, 1], [-5, -13, 9], [5, 15, -5], [15, 5, -5]}, {[-15, -5, -5], [-15, 5, 5], [15, -7, -1], [15, 7, 1]}, {[-15, -5, -5], [-5, -15, -5], [5, 13, 9], [15, 7, 1]}, {[-15, -5, 5], [-15, 5, -5], [15, -7, 1], [15, 7, -1]}, {[-15, -5, 5], [-5, -15, 5], [5, 13, -9], [15, 7, -1]}, {[-15, -3, -3], [-15, 3, 3], [15, -3, 3], [15, 3, -3]}, {[-15, -3, 3], [-15, 3, -3], [15, -3, -3], [15, 3, 3]}, {[-15, 5, -5], [-5, 15, -5], [5, -13, 9], [15, -7, 1]}, {[-15, 5, 5], [-5, 15, 5], [5, -13, -9], [15, -7, -1]}, {[-15, 7, -1], [-5, 13, -9], [5, -15, 5], [15, -5, 5]}, {[-15, 7, 1], [-5, 13, 9], [5, -15, -5], [15, -5, -5]}, {[-13, -6, -1], [0, -14, -11], [0, 11, 4], [13, 4, 5]}, {[-13, -6, 1], [0, -14, 11], [0, 11, -4], [13, 4, -5]}, {[-13, -4, -5], [0, -11, -4], [0, 14, 11], [13, 6, 1]}, {[-13, -4, 5], [0, -11, 4], [0, 14, -11], [13, 6, -1]}, {[-13, 4, -5], [0, -14, 11], [0, 11, -4], [13, -6, 1]}, {[-13, 4, 5], [0, -14, -11], [0, 11, 4], [13, -6, -1]}, {[-13, 6, -1], [0, -11, 4], [0, 14, -11], [13, -4, 5]}, {[-13, 6, 1], [0, -11, -4], [0, 14, 11], [13, -4, -5]}}

The radii above are:
2.100420126, 2.102918181, 2.104995276, 2.106740650, 2.108221156, 2.108878475, 2.109487662

I don't post the code now because it is not clean & tested.

Increasing Digits (e.g. to 25) will solve the problem.
(There are also options such as 'abserr', 'relerr').
BTW, why are you so sure that rkf45 is the best method for your ODE? If you are not, let dsolve choose the method.

Download Newton2.mw

For functions of a single variable such as fsolve(f(x), x=a..b)  the initial iteration is usually indeed (a+b)/2.
It is not easy to grasp the fsolve code for such a complex command. But you can do some tests to see the evaluation points of the function.
For example:

f:=proc(x)  print(x) ; cos(x)-x end:
fsolve(f, 0.5 .. 10);
                          5.250000000
                          5.250005250
                          5.249947500
                          2.875000000
                          2.875002875
                          2.874971250
                          7.625000000
                          7.625007625
                          7.624923750
                        3.87708345542864
                          1.687500000
                          1.687501688
                          1.687483125
                       0.782569730439221
                       0.742447770776506
                       0.739097995069612
                       0.739085133328356
                  0.7390851333290950851333284
                  0.7390851332151606416581575
                          0.7390851332

The restart command resets the "seed" of the random number generator to a default value.
You can remove the restart  or  insert a  randomize();  after it.

By default the arrow width depends on length. So, use e.g.:

PlotVector( [...], width=0.1);

psol := dsolve([ode1, ode2, x(0) = 50, y(0) = 10], numeric, output=listprocedure);
px:=eval(x(t),psol);

fsolve(px(t)=60, t=0..5);
       0.9330351416

You can also use
fsolve(px-60, 0..5);

Do you want to use Maple as a typewriter? Maple is supposed to be used as a powerful mathematical tool. However:

`2+3+10`;
                             2+3+10
% = parse(%);
                          2+3+10 = 15
 

BesselJ  is a "nice" function, so, Int should not have problems.
If you want a better timing and do not need >15 digits, try

evalf( Int(S, t=0..1, method=_d01ajc) );  # S = your sum

If it's still slow, extract the terms of the sum in a list
L:=convert(S,list);

and integrate them one by one to find the problematic ones.

for i to nops(L) do
   'L'[i]=evalf( Int(L[i], t=0..1, method=_d01ajc) );
od;

For non-square matrices having polynomial entries in a single variable, the best tool is the Smith Form.
So, the rank can be obtained at once, and it's easy to include everything in a procedure.

A := <b,1,3 | 4,b,6>:
B := <b,1,3 | 4,b+1,6>:
with(LinearAlgebra):
SmithForm(A), SmithForm(B);

So A has rank < 2 iff b=2, and B has always rank 2.

I have posted a year ago a work with this title at the Maple Application Center here.
Please do not expect miracles. It works mainly for not very complicated rational functions f. It has been reported that a 64 bit version of Maple is needed.