In[1]:= f[x_, y_] := Cos[x + y]
We need rho = (2 + n over n) terms to describe multinomial in 2 variables of degree <= n
In[3]:= n = 1
Out[3]= 1
In[6]:= ρ = Binomial[n + 2, n]
Out[6]= 3
In[7]:= pointsA = {
{0, 0},
{0, π/4},
{π/4, 0}
};
In[33]:= pointsB = {
{0, 0},
{π/8, π/8},
{π/4, π/4}
};
Set up Vandermonde matrix
In[8]:= VanderMondeLine[point_] := Block[{x = point[[1]], y = point[[2]]},
Return[{1, x, y}]
]
In[35]:= VA = VanderMondeLine /@ pointsA
VB = VanderMondeLine /@ pointsB
Out[35]= {{1,0,0},{1,0,π/4},{1,π/4,0}}
Out[36]= {{1,0,0},{1,π/8,π/8},{1,π/4,π/4}}
In[14]:= a = {a0, ax, ay};
In[37]:= ffA = f @@@ pointsA
ffB = f @@@ pointsB
Out[37]= {1,1/Sqrt[2],1/Sqrt[2]}
Out[38]= {1,1/Sqrt[2],0}
In[39]:= Solve[VA.a == ffA, a]
intA = VanderMondeLine[{x, y}].a //. %[[1]]
Out[39]= {{a0->1,ax->(2 (-2+Sqrt[2]))/π,ay->(2 (-2+Sqrt[2]))/π}}
Out[40]= 1+(2 (-2+Sqrt[2]) x)/π+(2 (-2+Sqrt[2]) y)/π
In[44]:= VB.a
Out[44]= {a0,a0+(ax π)/8+(ay π)/8,a0+(ax π)/4+(ay π)/4}
In[43]:= Solve[VB.a == ffB, a]
Out[43]= {}
In[25]:= ContourPlot[
f[x, y],
{x, 0, π/4},
{y, 0, π /4}
]

Out[25]= 
In[32]:= ContourPlot[
intA,
{x, 0, π/4},
{y, 0, π /4}
]
Out[32]=