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

Question

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}]

Answer 1

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]] *)