function [alpha,Psi_s,Psi_p,PsiSqNorm,P] = Hermite_PC(M,p_order)

n=10;
PRODOTTO2=zeros(n,n);
syms x1 x2 x3
a1=1.*x1.^0.*x2.^0.*x3.^0;
a2=x1;
a3=x2;
a4=x3;
a5=x1.^2-1;
a6=x2.^2-1;
a7=x3.^2-1;
a8=x1.*x2;
a9=x1.*x3;
a10=x2.*x3;
valori=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10];
for q=1:n
    for w=1:n
        h=matlabFunction(valori(q).*valori(w).*exp(-0.5.*x1.^2).*exp(-0.5.*x2.^2).*exp(-0.5.*x3.^2));
        PRODOTTO2(q,w)=round((1/(sqrt((2*pi)^3))*triplequad(h,-5,5,-5,5,-5,5,1.e-1)));
    end
end

n=10;
PRODOTTO3=zeros(n,n,n);
syms x1 x2 x3
a1=1.*x1.^0.*x2.^0.*x3.^0;
a2=x1;
a3=x2;
a4=x3;
a5=x1.^2-1;
a6=x2.^2-1;
a7=x3.^2-1;
a8=x1.*x2;
a9=x1.*x3;
a10=x2.*x3;
valori=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10];
for q=1:n
    for w=1:n
        for e=1:n
        h=matlabFunction(valori(q).*valori(w).*valori(e).*exp(-0.5.*x1.^2).*exp(-0.5.*x2.^2).*exp(-0.5.*x3.^2));
        PRODOTTO3(q,w,e)=round((1/(sqrt((2*pi)^3))*triplequad(h,-5,5,-5,5,-5,5,1.e-1)));
        end
    end
end