Bivariate(2D) Chebyshev implementation need to be checked because it doesn't work

Hello,

I have to implement the 2D-Chebyshev interpolation for a random function on a random rectangle.

The Interpolation polynomial is :

with the coeff

where the '' means that the first and the last summand are halved.

x_k and y_k are the chebyshev nodes.

T_i is the first degree chebyshev polynomial.

I don't know why code doesn't work.

Could someone help me. Here is my code