Symbolic matrice or vectors in Mathematica?

Question

Suppose I would like to define an expression, containing and array or matrix of symbolic elements.

Like this

p[x_, y_] := Sum[Sum[a[[i, j]]*x^i*x^j, {j, 1, 3}], {i, 1, 3}]

But if I try to evaluate this function, I see, that indexing is not symbolic:

Not only it swearing, but also returns incorrect result (zero).

So, is it possible to force Mathematica to know, that a is symbolic matrix of 3x3 containing no defined values. I.e. a[[1,2]] should evaluate to itself.

UPDATE

Some more detailed example

In[55]:= (* our system of equations with single solution *)
y == x*10 + 1 && y == -x*5 + 12

Out[55]= y == 1 + 10 x && y == 12 - 5 x

In[56]:= Clear[x, y]

In[57]:= (* solving a system *)
sol = Solve[y == x*10 + 1 && y == -x*5 + 12, {x, y}]

Out[57]= {{x -> 11/15, y -> 25/3}}

In[58]:= (* assigning found solutions *)
x = (x /. sol)[[1]]
y = (y /. sol)[[1]]

Out[58]= 11/15

Out[59]= 25/3

Now suppose we want to compose x and y into single vector
Using "functional" approach

In[60]:= Clear[x, y]

In[61]:= sol = Solve[x[2] == x[1]*10 + 1 && x[2] == -x[1]*5 + 12, {x[1], x[2]}]

Out[61]= {{x[1] -> 11/15, x[2] -> 25/3}}

In[62]:= (* assigning found solutions *)
x[1] = (x[1] /. sol)[[1]]
x[2] = (x[2] /. sol)[[1]]

Out[62]= 11/15

Out[63]= 25/3

In[9]:= x[1] /. sol

Out[9]= {1/2, 1/2}

In[10]:= x[2] /. sol

Out[10]= {-(Sqrt[3]/2), Sqrt[3]/2}

Unfortunately, x is not a vector now

In[65]:= (* we can't type it *)
x

Out[65]= x

In[66]:= (* we can't calculate it's length which should be 2 *)
Length[x]

Out[66]= 0

Now will try to use list approach

In[68]:= Clear[x]

(* and we fail just in the beginning *)
sol = Solve[x[[2]] == x[[1]]*10 + 1 && x[[2]] == -x[[1]]*5 + 12, x]

During evaluation of In[71]:= Part::partd: Part specification x[[2]] is longer than depth of object. >>

During evaluation of In[71]:= Part::partd: Part specification x[[1]] is longer than depth of object. >>

During evaluation of In[71]:= Part::partd: Part specification x[[2]] is longer than depth of object. >>

During evaluation of In[71]:= General::stop: Further output of Part::partd will be suppressed during this calculation. >>

During evaluation of In[71]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

Answer 1

Here is one way to keep everything symbolic:

p[x_, y_] := Total@Flatten@Table[a[i, j] x^i y^j, {j, 1, 3}, {i, 1, 3}]

p[3, 4]

12 a[1, 1] + 48 a[1, 2] + 192 a[1, 3] + 36 a[2, 1] + 144 a[2, 2] + 576 a[2, 3] + 108 a[3, 1] + 432 a[3, 2] + 1728 a[3, 3]

Of course, you can assign values to the a[i,j] in order to get numerical values for the polynomial.

Answer 2

Will this do?

A = Table[a[i,j],{i,3},{j,3}];

Then you can have

p[x_, y_] := Sum[Sum[A[[i, j]]*x^i*x^j, {j, 1, 3}], {i, 1, 3}];

and you will have e.g.:

Dimensions[A]

{3,3}

This is essentially what bill s proposes, but I am wrapping the symbolic a[i,j] in a matrix A.