Jaime,
You did not indicate that your x and y vectors were the full length of the temperature vector. In that case the formatting is even easier:
In[1]:=
xvectl = {1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4,
5};
yvectl = {1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
In[2]:= data = MapThread[{{#1, #2}, #3} &, {xvectl, yvectl, ts}]
Out[2]= {{{1, 1}, 111}, {{2, 1}, 122}, {{3, 1}, 133}, {{4, 1},
144}, {{5, 2}, 211}, {{1, 2}, 222}, {{2, 2}, 233}, {{3, 2},
244}, {{4, 3}, 311}, {{5, 3}, 322}, {{1, 3}, 333}, {{2, 3},
344}, {{3, 4}, 411}, {{4, 4}, 422}, {{5, 4}, 433}, {{1, 4},
444}, {{2, 5}, 511}, {{3, 5}, 522}, {{4, 5}, 533}, {{5, 5}, 544}}
However, Daniel is correct, Since the polynomial is guaranteed to go through every point, your interpolation will almost certainly be useless if MMA can even generate any answer at all.
My suggestion is to use NonlinearModelFit (or its related function, FindFit). In that case you can specify a function of any order and find the best fit. You can graph your data and the fitted function together and iterate until you are happy with the fitted function. Now you will have a lower order "best fit" to the data. In that case you need to change my code above because the fit functions take an {x,y,z} list,
data = Transpose[{xvectl, yvectl, ts}]
{{1, 1, 111}, {2, 1, 122}, {3, 1, 133}, {4, 1, 144}, {5, 2, 211}, {1,
2, 222}, {2, 2, 233}, {3, 2, 244}, {4, 3, 311}, {5, 3, 322}, {1, 3,
333}, {2, 3, 344}, {3, 4, 411}, {4, 4, 422}, {5, 4, 433}, {1, 4,
444}, {2, 5, 511}, {3, 5, 522}, {4, 5, 533}, {5, 5, 544}}
You can also experiment with FindFormula but I am not sure how it behaves with large amounts of data. I hope this helps.
Regards,
Neil