Off[General::"spell1"]

A variable shift

This is just an example of Mathematica code for a matrix of a linear transformation in a space of polynomials. Here the transformation is a variable shift. Mathematical definition is   (T p) (x) = p(
          x + c) .  The command below will produce the matrix of this linear transformation relative to the standard basis for polynomials.

MatrixForm[
      mS = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[
(IdentityMatrix[
                              5][[
                            k]] . {1, x, x^2, x^3, 
                            x^4}) /. {x -> (x + 
                              c)}
], x], x], 5], {k, 1, 
              5}]]]

Let me explain the above command, starting from inside.

I like to think of polynomials as "dot products" of a list of coefficients and powers:

{a0, a1, a2, a3} . {1, x, x^2, x^3}

This is the reason why I write the powers of the independent variable in such a strange way:

(IdentityMatrix[
            5][[3]] . {1, x, x^2, x^3, 
          x^4})

Above,  IdentityMatrix[5] is the 5-by-5 identity matrix, ,  IdentityMatrix[
              5][[3]] is its third row.

IdentityMatrix[5]

IdentityMatrix[
        5][[3]]

The following part makes the substitution

(IdentityMatrix[
              5][[3]] . {1, x, x^2, 
            x^3, x^4}) /. {x -> (x + c)}

{a0, a1, a2, a3} . {1, x, x^2, 
          x^3} /. {x -> (((x + c))^2)}

Now try

Expand[{a0, a1, a2, a3} . {1, x, x^2, 
            x^3} /. {x -> (((x + c))^2)}]

Now try, notice that I have to indicate the variable in the command below:

Collect[
      Expand[{a0, a1, a2, a3} . {1, x, x^2, 
              x^3} /. {x -> (((x + c))^2)}], x]

Now try, notice that I have to indicate the variable in the command below:

CoefficientList[
      Collect[Expand[{a0, a1, a2, a3} . {1, x, x^2, 
                x^3} /. {x -> (((x + c))^2)}], x], x]

The next command makes sure that I have all the coefficients. Notice the difference

CoefficientList[x + x^2, x]

PadRight[CoefficientList[x + x^2, x], 5]

Notice the index  k  in IdentityMatrix[5][[k]]; this index selects the k-th row. Since this matrix has 5 rows I go through all of them by using     Table[ -stuff here- ,{k,1,5}]

As often in Mathematica I obtain rows what should be columns, so I transpose.  I name the matrix  mS  and then show it easy to read using  MatrixForm[ ]

This is the end of the explanation.

The above process can be automated for any degree:

deg = 7

To create the powers of  x  I use table, but I have to prepend 1

Prepend[Table[(x^j), {j, 1, deg}], 1]

This is the only change in the above code.

deg = 9;  
    MatrixForm[
      mS = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[
(IdentityMatrix[
                              deg + 1][[
                            k]] . 
                          Prepend[Table[(x^j), {j, 1, deg}], 
                            1] ) /. {x -> (x + 
                              c)}
], x], x], 
              deg + 1], {k, 1, deg + 1}]]]

Legendre polynomials

Compared to the previous example, the only difference in the formula below is the definition of the transformation. In Mathematica we take derivative like this

D[x^3 + x^5, x]

and the second derivative is just

D[x^3 + x^5, {x, 2}]

Here is the differential operator

((1 - (x^2)) D[#, {x, 2}] - 2 x D[#, x]) &[
      (x^3)]

or,

((1 - (x^2)) D[#, {x, 2}] - 2 x D[#, x]) &[
      a0 + a1 x + a2 x^2]

Notice that in the above command a polynomial comes at the end in square brackets; 

& sign is just a part of the code; 

 the operator is defined inside of   ( -differential operator - )&  
 
 the symbol   #  is replaced by the polynomial given at the end in [  ]

MatrixForm[
      mL1 = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((1 - 
                                    (x^2)) D[#, {x, 2}] - 
                              2 x D[#, 
                                  x]) &[
(IdentityMatrix[5][[
                              k]] . {1, x, x^2, x^3, 
                              x^4})
]
]], x], x], 5], {k, 1, 5}]]]

esL1 = Eigensystem[mL1]

Reverse[({1, x, x^2, x^3, x^4} . #) & /@ 
        esL1[[2]]]

Normalization is by p(1) = 1.

\(\(\((#\/\(# /. {x -> 1}\))\) &\) /@ 
      Reverse[\(\(({1, x, x^2, x^3, x^4} . #)\) &\) /@ 
          esL1[\([2]\)]]\)

Plot[{1, x, (1/2) (-1 + 3 (x^2)), 
        (1/2) (-3 x + 5 (x^3)), 
        (1/8) (3 - 30 (x^2) + 35 (x^4))}, {x, -1.1, 1.1}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}}, 
      PlotRange -> {-1.2, 1.2}, ImageSize -> 400]

Expand[
      Table[(1/((2^n) n!)) D[(((x^2 - 1))^n), {x, n}], {n, 0, 
          4}]]

Clear[ps]; 
ps = 
      Expand[Table[(1/((2^n) n!)) D[(((x^2 - 1))^n), {x, n}], {n, 
            0, 8}]]

Table[
        Integrate[
          ps[[j]] ps[[
              k]], {x, -1, 1}], {j, 1, 
          Length[ps]}, {k, 1, Length[ps]}] // MatrixForm

Legendre polynomials (set your degree)

deg = 7;

Clear[mL1]; 
    MatrixForm[
      mL1 = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((1 - 
                                    (x^2)) D[#, {x, 2}] - 
                              2 x D[#, 
                                  x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esL1 = Eigensystem[mL1]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esL1[[2]]]

Normalization: p(1) = 1

\(polysL1 = \(\((#\/\(# /. {x -> 1}\))\) &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esL1[\([2]\)]]\)

Plot[Evaluate[polysL1], {x, -1.1, 1.1}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}, {Yellow, 
            Thickness[0.007]}}, PlotRange -> {-1.5, 1.5}, 
      ImageSize -> 400]

Expand[
      Table[(1/((2^n) n!)) D[(((x^2 - 1))^n), {x, n}], {n, 0, 
          deg}]]

Laguerre polynomials

deg = 7;  
    MatrixForm[
      mL2 = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
(x D[#, {x, 
                                    2}] + (1 - x) 
                                D[#, x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esL2 = Eigensystem[mL2]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esL2[[2]]]

Noramlization: p(0) = 1.

\(polysL2 = \(\((#\/\(# /. {x -> \ 0}\))\) &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esL2[\([2]\)]]\)

Plot[Evaluate[polysL2], {x, -1, 5}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}}, 
      PlotRange -> {-2, 4}, ImageSize -> 400]

Bessel polynomials

deg = 7;  
    MatrixForm[
      mB = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((x^2) D[#, {x, 
                                    2}] + 
                              2 (1 + x) 
                                D[#, x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esB = Eigensystem[mB]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esB[[2]]]

Normalization: p(0) = 1

\(polysB = \(\((#\/\(# /. {x -> 0}\))\) &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esB[\([2]\)]]\)

Plot[Evaluate[polysB], {x, -1, 1}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}}, 
      PlotRange -> {-2, 5}, ImageSize -> 400]

Chebyshev polynomials

deg = 5;

Clear[mCh]; 
    MatrixForm[
      mCh = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((1 - 
                                    (x^2)) D[#, {x, 2}] - (x) 
                                D[#, x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esCh = Eigensystem[mCh]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esCh[[2]]]

Normalization: p(1) = 1

\(polysCh = \(\((#\/\(# /. {x -> 1}\))\) &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esCh[\([2]\)]]\)

Plot[Evaluate[polysCh], {x, -1.1, 1.1}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}, {Yellow, 
            Thickness[0.007]}}, PlotRange -> {-1.5, 1.5}, 
      ImageSize -> 400]

Table[TrigExpand[Cos[n ArcCos[x]]], {n, 0, 5}]

Chebyshev polynomials - many

deg = 10;

Clear[mCh]; 
    MatrixForm[
      mCh = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((1 - 
                                    (x^2)) D[#, {x, 2}] - (x) 
                                D[#, x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esCh = Eigensystem[mCh]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esCh[[2]]]

Normalization: p(1) = 1

\(polysCh = \(\((#\/\(# /. {x -> 1}\))\) &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esCh[\([2]\)]]\)

Plot[Evaluate[polysCh], {x, -1.1, 1.1}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}, {Yellow, 
            Thickness[0.007]}}, PlotRange -> {-1.5, 1.5}, 
      ImageSize -> 400]

Table[TrigExpand[Cos[n ArcCos[x]]], {n, 0, deg}]

Looking for a general formula for Chebyshev polynomials

Even degree

nn = 5;  
    ttt = Table[{(((2 nn))^2) - (((2 k))^2),  (2 k + 2) (2  k + 
                1)}, {k, 0, nn - 1}]

\(\(\((#[[2]]\/#[[1]])\) &\) /@ ttt\)

(5/2) 512

(7/8) 1280

2^(2  nn - 1) 
      Table[(((-1))^(j + 1)) 
          Product[ (((2 k + 2) (2  k + 1))/((((2 nn))^2) - (((2 k))^2))), {k, j, nn - 1}], {j, 0, nn}]

nn = 11;

2^(2  nn - 1) 
      Table[(((-1))^(j + 1)) 
          Product[ (((2 k + 2) (2  k + 1))/((((2 nn))^2) - (((2 k))^2))), {k, j, nn - 1}], {j, 0, nn}]

Total[%]

Hermite polynomials

deg = 5; 
    MatrixForm[
      mH = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((1) D[#, {x, 
                                    2}] - (2 x) 
                                D[#, x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esH = Eigensystem[mH]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esH[[2]]]

The normalization: the leading coefficient is (2^n).

\(polysH = \(Expand[\((\(2\^\(Length[CoefficientList[#, x]] - 
                      1\)\) #\/Last[CoefficientList[#, x]])\)] &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esH[[2]]]\)

Plot[Evaluate[polysH], {x, -3, 3}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}, {Yellow, 
            Thickness[0.007]}}, PlotRange -> {-35, 40}, 
      ImageSize -> 400]

Expand[
      Table[Simplify[(((-1))^n) Exp[x^2] 
            D[Exp[-x^2], {x, n}]], {n, 0, deg}]]

Experiment

deg = 3;

Clear[mE]; 
    MatrixForm[
      mE = Transpose[
          Table[PadRight[
              CoefficientList[
                Collect[
                  Expand[Simplify[
((1 + 
                                    0  x) D[#, {x, 2}] + (1 + x) D[#, 
                                  x]) &[
(IdentityMatrix[deg + 1][[
                              k]] . 
                            Prepend[Table[(x^j), {j, 1, deg}], 
                              1])
]
]], 
                  x], x], deg + 1], {k, 1, deg + 1}]]]

esE = Eigensystem[mE]

Reverse[(Prepend[Table[(x^j), {j, 1, deg}], 1] . #) & /@ 
        esE[[2]]]

Normalization: p(1) = 1

\(polysE = \(\((#\/\(# /. {x -> 1}\))\) &\) /@ 
        Reverse[\(\((Prepend[Table[x\^j, {j, 1, deg}], 1] . #)\) &\) /@ 
            esE[\([2]\)]]\)

Plot[Evaluate[polysCh], {x, -1.1, 1.1}, 
      PlotStyle -> {{Red, Thickness[0.007]}, {Cyan, 
            Thickness[0.007]}, {Green, Thickness[0.007]}, {Magenta, 
            Thickness[0.007]}, {Blue, Thickness[0.007]}, {Yellow, 
            Thickness[0.007]}}, PlotRange -> {-1.5, 1.5}, 
      ImageSize -> 400]

Table[TrigExpand[Cos[n ArcCos[x]]], {n, 0, 5}]