How to calculate the mean and the variance for the function with Multidimensional Polynomial Chaos Approach?

Question

I would like to calculate the mean and the variance for the function

f(x1,x2)=pi*(x1-1)*sin(pi*x1)*(1-x2^2); (x1,x2)->Uniform distribution [-1,1]

using the multidimensional Polynomial Chaos Approach.

First thing that I did is to normalize the Legendre Polynomials as Pi(x)/sqrt(1/(2*i+1)).

Second thing, define the points and weights from the Gauss-Legendre quadrature for m+1.

Third thing, make the tensor product or all combination of the points and weights.

...and solve the system for the unknown coefficients.

I take the degree p=4, P=(4+2)!/(4!*2!)=15 unknown coefficients, full quadrature scheme 25 calculation. I need to get the result for the mean 0.66649 and for the variance 2.7434541 but somewhere somehow something is wrong and I get the mean 0.666 and the variance 1.05 instead 2.74344541.

Any help, comments? Thanks.

Matlab script

%% Order of Polynomial p=4;

%% Number of nodes m=p+1;

%% Calls the gauss-legendre quadrature [x,w]=lgwt(m,-1,1); xq=x; wq=w/2;

%% number of unknows d=2;

%% unknown coefficients P=factorial(p+d)/(factorial(p)*factorial(d));

%% zeros of function yuq=zeros(P+1,1);

%% All combination comx=allcomb(xq,xq); comw=allcomb(wq,wq);

[r,c]=size(comx);

a=comx(:,1:round(c/2)); b=comx(:,round(c/2)+1:end);

s=length(comx);

%% Non Intrusive gPC Solver

for i=0:P

for j=1:s

    a1=a(j,1);
    b1=b(j,1);

    yuq(i+1,1)=yuq(i+1,1)+pi*(a1-1)*sin(pi*a1)*(1-b1^2)*prod(legendreP2(i,comx(j,:)))*prod(comw(j,:));

end

end

%% Mean and Varinace calculator

m=yuq(1); v=sum(yuq(2:end,1).^2);