Pdf1= 2 Sqrt[2/[Pi]] (x1-x2) exp(1/2 (-x1 (3 x1-x2)-x2 (3 x2-x1)))

rg= {{x1, -Infinity, Infinity}, {x2, -Infinity, x1}};
ContourPlot[Pdf1, Sequence @@ rgn // Evaluate,  ImageSize -> Small]

Clear[int]; int[a_, b_] :=
 Integrate[ a b  Pdf1, Sequence @@ rg // Evaluate]

Clear[nint]; nint[a_, b_] :=
 NIntegrate[ a b  Pdf1, Sequence @@ rg // Evaluate]

p = 2; pol0 =
 Table[Table[x1^i x2^(p1 - i), {i, 0, p1}] // Flatten // Union,
    {p1, 0, p}] // Flatten // Join

pol = Orthogonalize[pol0, int[#1, #2] &]

Map[ContourPlot[# Pdf1  , Sequence @@ rgn // Evaluate,
   PlotPoints -> 15, PlotRange -> All] &, pol]

pol = Orthogonalize[pol0, intn[#1, #2] &,  Method -> "Reorthogonalization"];

p = 6; pol0 =
 Table[Table[x1^i x2^(p1 - i), {i, 0, p1}] // Flatten // Union,
    {p1, 0, p}] // Flatten // Join;
pol = Orthogonalize[pol0, intn[#1, #2] &,  Method -> "Reorthogonalization"];

gs[vecs_, ip___] := Module[{ovecs = vecs},
Do[ovecs[[i]] -= Projection[ovecs[[i]], ovecs[[j]], ip], {i, 2,
Length[vecs]}, {j, 1, i - 1}]; ovecs];
pol1 = gs[pol0, Function[
    NIntegrate[ ##  Pdf1, Sequence @@ rg // Evaluate]]];

mat = ParallelTable[int[i, j], {i, pol0}, {j, pol0}];
   eigs = Eigensystem[mat // N];

Map[ContourPlot[#  Pdf1, Sequence @@ rgn // Evaluate,
   PlotPoints -> 15, PlotRange -> All] &,
 pol = eigs[[2]].pol0/Sqrt[eigs[[1]]] // Chop]