NN  =  4;

  

  

Clear[xi, 
      xi0]
N[
        Solve[D[LegendreP[NN, x], x] == 
            0]];
xi0 = 
        Sort[Union[Re[x /. %], {1, -1}], 
          Smaller];
Do[
        xi[i] = xi0[[i + 1]], {i, 0, NN}];
xi0

Clear[
      wi]
Do[
        wi[i] = (2/(NN (NN + 1))) (1/LegendreP)[NN, xi[i]]^2, {i, 0, 
          NN}];
wi0 = Table[wi[i], {i, 0, NN}]

Do[
        M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[M[n, n] = 1, {n, 0, NN}];
matMid = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMid]

  
Do[M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[
        M[n, n - 1] = (n/(2  n - 1)), {n, 1, 
          NN}];
Do[
        M[n, n + 1] = ((n + 1)/(2  n + 3)), {n, 1, 
          NN}];
matMml = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMml]

  
Do[M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[
        M[n, p] = 2  n + 1, {n, 0, NN}, {p, n + 1, NN, 
          2}];
matMD1 = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMD1]

  
Do[M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[
        M[n, p] = (n + (1/2)) (p (p + 1) - n (n + 1)), {n, 0, 
          NN}, {p, n + 2, NN, 2}];
matMD2 = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMD2]

(* 
      Operator (((d^2) f)/(dx^2))  - 4  (df/dx)  + 
        4  f  
   matL  :  keeps as a List (matrix)  
L[i, j] :  
              Keeps as components       *) 

matL  =  
        matMD2  -  4 matMD1  +  4  matMid;
Do[
        L[i, j] = matL[[i + 1, j + 1]], {i, 0, NN}, {j, 0, 
          NN}];
MatrixForm[matL]
Eigenvalues[matL]

Do[
        fm[n] = LegendreP[n, -1], {n, 0, 
          NN}];
Do[fp[n] = LegendreP[n, 1], {n, 0, NN}];
Do[L[NN - 1, n] = fm[n], {n, 0, NN}];
Do[L[NN, n] = fp[n], {n, 0, NN}];
matL = 
        Table[L[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matL]

Clear[
      fsol]
c0 =  -4 (Exp[x]  -  Sinh[1]/Sinh[2] Exp[2  x]) /. 
          x -> -1;
fsol[x_]  =  
        Exp[x]  -  Sinh[1]/Sinh[2] Exp[2  x] + 
          (c0/4);
N[fsol[1]]

(* R . H . S .   
n = 
        NN - 1 and NN lines are for boundary condition of Tau method  *) 
 (*  Make spectral coeffcients of RHS*) 

Do[
      rhs[i] = N[Exp[xi[i]] + c0], {i, 0, NN}]

Do[
        gam[n] = Sum[LegendreP[n, xi[i]]^2 wi[i], {i, 0, NN}], {n, 0, 
          NN}];
Do[
        coefrhs[n] = (1/gam)[n] 
            Sum[rhs[i] LegendreP[n, xi[i]] wi[i], {i, 0, NN}], {n, 0, 
          NN}];
coefrhs[NN - 1] = 0;
coefrhs[NN] = 0;
vecRHS  =  Table[coefrhs[n], {n, 0, NN}]

(*  Solve for the spectral coefficients  *) 
matLN  =  N[matL];
vecSOL  =  
        LinearSolve[matLN,  vecRHS];
Do[
        coeff[n] = vecSOL[[n + 1]], {n, 0, 
          NN}];

(* 
      Compute Interpolant of f (x)  *) 

interpf[x_]  = 
      Sum[coeff[n] LegendreP[n, x], {n, 0, NN}]

(* 
      Plot interpolant and exact solution  *) 

Plot[{interpf[x], fsol[x]}, {x, -1, 
        1}]

(*  USE CHEBYSHEV !!  *)

NN  =  
        4;
  

  

Clear[xi, xi0]
Do[
        xi[i] = Cos[((π i)/NN)], {i, 0, 
          NN}];
xi0 = Table[xi[i], {i, 0, NN}]

Clear[
      wi]
Do[wi[i] = (π/NN), {i, 1, NN - 1}];
wi[0] = (π/(2  NN)); 
    wi[NN] = (π/(2  NN));
wi0 = Table[wi[i], {i, 0, NN}]

Do[
        M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[M[n, n] = 1, {n, 0, NN}];
matMid = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMid]


      
Do[M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[
        M[n, n - 1] = (1/2) (1 + KroneckerDelta[0, n - 1]), {n, 1, 
          NN}];
Do[M[n, n + 1] = (1/2), {n, 1, NN}];
matMml = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMml]


      
Do[M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[
        M[n, p] = (2/((1 + KroneckerDelta[0, n]))) p, {n, 0, NN}, {p, 
          n + 1, NN, 2}];
matMD1 = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMD1]


      
Do[M[i, j] = 0, {i, 0, NN}, {j, 0, NN}];
Do[
        M[n, p] = (1/((1 + KroneckerDelta[0, n]))) 
            p ((p^2) - (n^2)), {n, 0, NN}, {p, n + 2, NN, 
          2}];
matMD2 = 
        Table[M[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matMD2]

(* 
      Operator (((d^2) f)/(dx^2))  - 4  (df/dx)  + 
        4  f  
   matL  :  keeps as a List (matrix)  

            
            L[i, j] :  
              Keeps as components       *) 

matL  =  
        matMD2  -  4 matMD1  +  4  matMid;
Do[
        L[i, j] = matL[[i + 1, j + 1]], {i, 0, NN}, {j, 0, 
          NN}];
MatrixForm[matL]
Eigenvalues[matL]

Do[
        fm[n] = LegendreP[n, -1], {n, 0, 
          NN}];
Do[fp[n] = LegendreP[n, 1], {n, 0, NN}];
Do[L[NN - 1, n] = fm[n], {n, 0, NN}];
Do[L[NN, n] = fp[n], {n, 0, NN}];
matL = 
        Table[L[i, j], {i, 0, NN}, {j, 0, NN}];
MatrixForm[matL]

Clear[
      fsol]
c0 =  -4 (Exp[x]  -  Sinh[1]/Sinh[2] Exp[2  x]) /. 
          x -> -1;
fsol[x_]  =  
        Exp[x]  -  Sinh[1]/Sinh[2] Exp[2  x] + 
          (c0/4);
N[fsol[1]]

(* R . H . S .   
n = 
        NN - 1 and NN lines are for boundary condition of Tau method  *) 
 (*  Make spectral coeffcients of RHS*) 

Do[
      rhs[i] = N[Exp[xi[i]] + c0], {i, 0, NN}]

    
Do[
        gam[n] = Sum[ChebyshevT[n, xi[i]]^2 wi[i], {i, 0, NN}], {n, 0, 
          NN}];
Do[
        coefrhs[n] = (1/gam)[n] 
            Sum[rhs[i] ChebyshevT[n, xi[i]] wi[i], {i, 0, NN}], {n, 0, 
          NN}];
coefrhs[NN - 1] = 0;
coefrhs[NN] = 0;
vecRHS  =  Table[coefrhs[n], {n, 0, NN}]

(*  Solve for the spectral coefficients  *) 
matLN  =  N[matL]
vecSOL  =  LinearSolve[matLN,  vecRHS]
Do[
        coeff[n] = vecSOL[[n + 1]], {n, 0, 
          NN}];

(* 
      Compute Interpolant of f (x)  *) 

interpf[x_]  = 
      Sum[coeff[n] ChebyshevT[n, x], {n, 0, NN}]

(* 
      Plot interpolant and exact solution  *) 

Plot[{interpf[x], fsol[x]}, {x, -1, 
        1}]