F:=(x,y)->sin(x)*cos(y);

M:= Matrix(2,3, (i,j) -> [i,2*i +j]);

map('t ->F(t[1], t[2])',M); 
evalf(%);

          [sin(1) cos(3)    sin(1) cos(4)    sin(1) cos(5)]
          [                                               ]
          [sin(2) cos(5)    sin(2) cos(6)    sin(2) cos(7)]

Your function a - please let us call it F - is defined on 'pairs', i.e. in IR^2.
So you have to provide that as input.

  F:=(x,y)->sin(x)*cos(y);

  [seq( [i,i/2], i = 0 .. 360, 45)];
  map('t ->F(t[1], t[2])',%);

I am a bit slow with typing, but as I already did here is my answer 
(even if now partially repeating others ...)

The help page for 'sum' recommends to use 'add':

  add(f(t), t = -2 .. 2);
                                  6

If you want to use 'sum' you can use quotation marks:

  sum('f'(t), t = -2 .. 2);
                                  6


Note however if using the upper case, inert, version it will fall
back to your initial problem if you are going to evaluate it:

  Sum('f'(t), t = -2 .. 2); value(%);

    Error, (in f) cannot determine if this expression is ...


You have to use 2 quotation marks then:

  Sum(''f''(t), t = -2 .. 2); value(%);

                                  6

A better way is to use some 'piecewise' definition:

  F:= i -> piecewise( 0<i, i, -i);
  sum(F(t), t = -2 .. 2);

                                  6

Compare the different behaviours:

  f(t);

    Error, (in f) cannot determine if this expression is ...

  F(t);

                        { t           0 < t
                        {
                        { -t        otherwise


The latter is possibly what you intend for your original F.


In any case it would be good to care for 'undecidable', as the
system does not work like it would do in common programming
(but do it with care for intended use):

  G:= i -> piecewise( 0<i, i, i<=0, -i, `doesNotWork`(i));

It does not know a priori that you want numerical inputs only
and would accept other stuff as well (where ' t < 0 ' makes no
sense, like complex numbers, matrices or functions) and 'sum'
does just that, contrary to 'add'.


Also you can provide the 'data type' like

  g:= proc(i::integer)
    if 0 < i then return i else return -i end if 
  end proc;

Using it will show the difference between sum and add.


A final example regarding the use of F and G from above using
complex numbers and variables:

  1/(x^2+1); convert(%, fullparfrac);

This gives you the partial fraction for 1/(x^2+1) in symbolic
form using the RootOf construct and it is a legal input for sum:

  'sum( F(-1/2*_alpha/(x-_alpha)), _alpha = RootOf(_Z^2+1))';
  allvalues(%); simplify(%);

                       - 1/(x^2+1)

Observe the minus sign: This is correct, since the inputs for F
are not Reals, especially can not be identified to be larger then
zero, so F takes the negative value.

Now do the same, but using G instead of F ... and after that do

  eval(%, `doesNotWork` =' z -> z');  simplify(%);

Up to a stupid simplification problem you get 1/(x^2+1).


Hope it helps: you have to care for intended use.

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