# How to show the Chebyshev polynomials in the solution of a differential equation?

In the following code I solved a differential equation on Mathematica 11. The solution has two arbitrary constants, which, if set to 0 and 1, in previous versions of Mathematica used to give as the solution the Chebyshev polynomials.

How can I tell to Mathematica that I would like to set the constants in order to keep (and later plot) the part of the solution that represents the Chebyshev polynomials?

sol4[n_] := y[x] /. DSolve[ (1 - x^2) y''[x] - x y'[x] + n^2 y[x] == 0, y[x], x] [[1]] /. {C[1] ->1, C[2]->0}

lista = Table[ sol4[i], {i, 1, 5}];
gra1 = Plot[lista, {x, -1, 1}]

listaorig = Table[ ChebyshevT[i,x], {i, 1, 5}];
gra2 = Plot[ listaorig, {x, -1, 1}]

• TIP: use DSolveValue instead of DSolve, i.e., DSolveValue[(1 - x^2) y''[x] - x y'[x] + n^2 y[x] == 0, y[x], x] /. {C[1] -> 1, C[2] -> 0} – AccidentalFourierTransform Apr 29 '18 at 13:27

Clear[sol4]

sol4[n_, x_] =
DSolveValue[(1 - x^2) y''[x] - x y'[x] + n^2 y[x] == 0, y[x],
x] /. {C[1] -> c1, C[2] -> c2} // FullSimplify[#, -1 < x < 1] &

(* c1 Cosh[n ArcTanh[x/Sqrt[-1 + x^2]]] + c2 Sin[n ArcSin[x]] *)


Solve for the constants such that the sol4[n, x] is equal to ChebyshevT[n,x]. Since there are two unknowns, you need a second equation: equate the derivatives.

Assuming[-1 < x < 1,
Solve[{sol4[n, x] == ChebyshevT[n, x],
D[sol4[n, x], x] == D[ChebyshevT[n, x], x]}, {c1, c2}][[1]] //
FullSimplify]

(* {c1 -> Cos[(n π)/2], c2 -> Sin[(n π)/2]} *)


After clearing the definition, redefine the function with these constants

Clear[sol4]

sol4[n_, x_] =
DSolveValue[(1 - x^2) y''[x] - x y'[x] + n^2 y[x] == 0, y[x],
x] /. {C[1] -> Cos[n π/2], C[2] -> Sin[n π/2]} //
FullSimplify[#, -1 < x < 1] &

(* Cos[n ArcCos[x]] *)


Verify that sol4[n, x] == ChebyshevT[n, x]

sol4[n, x] == ChebyshevT[n, x] // FullSimplify

(* True *)


The plots are as expected

lista = Table[sol4[n, x], {n, 1, 5}];
gra1 = Plot[lista, {x, -1, 1}]


EDIT: Alternatively, looking at the plots for ChebyshevT[n,x], just include the values for y[0] and y[1] to define the constants.

Clear[sol4]

sol4[n_, x_] =
DSolveValue[{(1 - x^2) y''[x] - x y'[x] + n^2 y[x] == 0,
y[0] == Cos[n π/2], y[1] == 1}, y[x], x] //
FullSimplify[#, -1 < x < 1] &

(* Cos[n ArcCos[x]] *)