Your Y(.01) is a list, you can access elements by their indices

These elements are equation and you can access the right and left hand side

a=rhs( Y(.01)[2] );

                       a = 0.989950166993563752

b=rhs( Y(.01)[3] );

                       b = -1.00994983384860371

If you want to define a and b as such you use double points
 

More or less you should find out what you need from the help
pages for dsolve/numeric (it needs some time to get used to
the styles).

For example:

  dsol2 := dsolve({ode,con}, numeric, output=listprocedure);

  dsol2(0.01);
  a,b:=op(2..3,%);

  a, b := y(x)(0.01) = 0.989950166993563752,

        /d      \
        |-- y(x)|(0.01) = -1.00994983384860371
        \dx     /

The last statement assigns for a sequence while the 'op' will
take the operants (= entries in list). Of course you can do it
separately for each entry as well.

Note that the values may depend on your setting of Digits.

The following procedure returns the transform matrix using rationals only
and works for larger n, say n=512. To compare with your original version
in tests however one would smaller values for n and comparing is better
done after using evalf on both matrices (to avoid symbolic troubles and
evalf is reasonable fast for Gamma, which will be used in the original
solution for larger n, 2 <= n):

  theA:=proc(n::posint)
  local A,i,j, q,r;

  A:=Array(0..n,0..n, fill=0);
  
  # first row
  A[0,0]:=1;
  for j from 2 by 2 to n do
    A[0,j]:= -1/(j^2-1);
  end do;
  
  # second row
  A[1,1]:=1;
  for j from 3 by 2 to n do
    A[1,j]:= -3/(j^2-4);
  end do;
  
  for i from 2 to n do
    # upper triangle: fill sub-diagonals
    # using recursion for F on those diagonals
    q:= (2*i+1)/(2*i-1);
    for j from i to n do
      if (i-j) mod 2 = 0 then
        r:= j/(j-1) * (i+j-2)/(i+j+1);
        A[i,j]:= A[i-1,j-1]*q*r; 
      end if;
    end do;
  end do;
  
  return convert(A,Matrix);
  end proc; #########################################################


For that I tried to use evalhf, but it needed endless time for larger
inputs. Since I seem to be too stupid to use compile for zero based
Arrays as input I re-coded it as Matrix version ( M[i,j]=A[i-1,j-1] )
for a floating point version, 2 <= n:

  theM_numerical:=proc(n::posint, M::Matrix(datatype = float[8]))
  local i,j, q,r;
  
  # first row
  M[1,1]:=1;
  for j from 3 by 2 to n+1 do
    M[1,j]:= -1.0/((j-1)*(j-1)-1.0);
  end do;
  
  # second row
  M[2,2]:=1;
  for j from 4 by 2 to n+1 do
    M[2,j]:= -3.0/((j-1)*(j-1)-4.0);
  end do;
  
  for i from 3 to n+1 do
    # upper triangle: fill sub-diagonals
    q:= (2.0*i-1.0)/(2.0*i-3.0);
    for j from i to n+1 do
      if (i-j) mod 2 = 0 then
        r:= (j-1.0)/(j-2.0);
        r:= r * (i+j-4.0)/(i+j-1.0); 
        M[i,j]:= M[i-1,j-1]*q*r; 
      end if;
    end do;
  end do;
  
  return NULL;
  end proc; #########################################################


This can be compiled (I needed some experiments to find out what that
wants as input type for the array):

  theM_cp:=Compiler:-Compile(theM_numerical);

And it works fast for n=512 (define the test value and provide the
matrix as argument for the compiled version):

  dTst:= 512;
  arrM:=Matrix(1..dTst+1,1..dTst+1, datatype = float[8] ):
  st:=time():
  theM_cp(dTst,arrM);
  `sec needed`=time() - st;

                          sec needed = 0.016

Now arrM is updated to hold the transform matrix (in double precision).

As speed is not all let us check for 'exactness' for Digits=14 (which
is appropriate here for me):

  'arrM - evalf(theA(dTst))'; 
  convert(%, listlist): ListTools[Flatten](%): convert(%,set); 
  fnormal(%); simplify(%,zero);

                              {0.}

So exactness seems to be satisfactory (I did not directly convert to
a set, just to see some intermediate results).

It is similar to acer's version (but needs no specific length and works for any input to give the constant function):

f:= i -> unapply(i,x);

or

g:= i -> (x -> i);

f(3);
                                x -> 3
f(3.5)(t);
                                 3.5
g(Pi)(exp(1));
                                  Pi

make it a polynomial expression by using convert(%, polynom)

or like this:

  F:= u -> u*D(u);

  F(y)(x);
                             y(x) D(y)(x)

  F(sin)(x);
                            sin(x) cos(x)

  G:= (u,v) -> sin@(u*v);

  G(u,v)(x);
                            sin(v(x) u(x))

  LegendreQ(1,x); convert(%, Int): combine(%);

  LegendreQ(1,x); convert(%,hypergeom); 
    convert(%, Int):  simplify(%); combine(%);

In both case you get the same integral representation

                         1
                        /          1/2
                       |        _t1
                       |   - ------------- d_t1
                       |          2
                      /      2 (-x  + _t1)
                        0

  value(%) assuming x::real;
  convert(%, piecewise,x);
 
                         -1 + x arctanh(1/x)

which you convert to a log or take the Taylor series

  series(%,x=0) assuming x::real;

                                     2        4      6
      - 1 + -1/2 I signum(x) Pi x + x  + 1/3 x  + O(x )

p3 = 0, p5 = 0, p7 = 0 and the classic interface responds by "Warning, solutions may have been lost"

I think it is a linear system and I would write it as such

fsolve(Int(Int(1,rbar=0..RootOf(f(r,z)-1,r)),z=0..1));

  a1:=Int(1/(2*A),x=-A..A);
  IntegrationTools[Change](a1,x=20*log[10](u)):
  combine(%); # <--- !
  op(1,%);

                                  10
                              ----------
                              A u ln(10)

  IntegrationTools[Change](a1,x=K*log[10](u), u )  assuming 1 < K;
  #                                          ^^^ provide new variable,
  #                                              you have 'u' and 'K'
  combine(%);
  op(1,%);

                                    K
                            1/2 ----------
                                ln(10) A u

test2 := proc(x::uneval)
local X;
X:=eval(x);
  x:=map('t -> t + 1', X);
  NULL;
end proc;

x:=10;
y:=20;
L:=[x,y];
test2(L);
x,y,L;

                           10, 20, [11, 21]

it would not accept your input

  2^3^4;

  Error, ambiguous use of `^`, please use parentheses

  (2^3)^4;
                                    4096

  2^(3^4);
                          2417851639229258349412352

So it was likeli interpreted as the latter case - bang

plot3d(-abs(x-y) , x=-1..1, y=-1..1
    #, gridstyle=triangular
    , axes=boxed
 );

You can increase Digits or use more numpoints as well, 
see ?plot3d/option

What you do is to write Euler's integral as 2F1, but it is not
continuous in its parameters (else it would be 1/2 in that case)
and Maple does not do it.

To prove it: 

  eq3; combine(%);
  op(1,%);
  limit(%,k=infinity);
                                  0

So the integrand is 0 in the limit. Now justify interchanging
limit and integral to get your desired result.

Not that complete as the package as described in the linked paper above (especially they do have various methods), but the following fits my basic needs (and I am too 'lazy' for a full solution, that's on Maple, if they want to compete ...):

Download 102_computing_Gauss_integration_rules.mws (30 kb)