> with(linalg);
Warning, new definition for norm
Warning, new definition for trace

  [BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp,

        Wronskian, addcol, addrow, adj, adjoint, angle, augment,

        backsub, band, basis, bezout, blockmatrix, charmat, charpoly,

        cholesky, col, coldim, colspace, colspan, companion, concat,

        cond, copyinto, crossprod, curl, definite, delcols, delrows,

        det, diag, diverge, dotprod, eigenvals, eigenvalues,

        eigenvectors, eigenvects, entermatrix, equal, exponential,

        extend, ffgausselim, fibonacci, forwardsub, frobenius,

        gausselim, gaussjord, geneqns, genmatrix, grad, hadamard,

        hermite, hessian, hilbert, htranspose, ihermite, indexfunc,

        innerprod, intbasis, inverse, ismith, issimilar, iszero,

        jacobian, jordan, kernel, laplacian, leastsqrs, linsolve,

        matadd, matrix, minor, minpoly, mulcol, mulrow, multiply,

        norm, normalize, nullspace, orthog, permanent, pivot,

        potential, randmatrix, randvector, rank, ratform, row, rowdim,

        rowspace, rowspan, rref, scalarmul, singularvals, smith,

        stackmatrix, submatrix, subvector, sumbasis, swapcol, swaprow,

        sylvester, toeplitz, trace, transpose, vandermonde, vecpotent,

        vectdim, vector, wronskian]

> A:=matrix(6,6,[-2,1,0,0,0,0,1,-2,1,0,0,0,0,1,-2,1,0,0,0,0,1,-2,1,0,0,0
> ,0,1,-2,1,0,0,0,0,1,-2]); # THIS IS THE MATRIX A

                    [-2     1     0     0     0     0]
                    [                                ]
                    [ 1    -2     1     0     0     0]
                    [                                ]
                    [ 0     1    -2     1     0     0]
               A := [                                ]
                    [ 0     0     1    -2     1     0]
                    [                                ]
                    [ 0     0     0     1    -2     1]
                    [                                ]
                    [ 0     0     0     0     1    -2]

> AINV:=inverse(A); # THIS IS THE INVERSE OF A

              [-6/7    -5/7     -4/7     -3/7     -2/7     -1/7]
              [                                                ]
              [-5/7    -10/7    -8/7     -6/7     -4/7     -2/7]
              [                                                ]
              [-4/7    -8/7     -12/7    -9/7     -6/7     -3/7]
      AINV := [                                                ]
              [-3/7    -6/7     -9/7     -12/7    -8/7     -4/7]
              [                                                ]
              [-2/7    -4/7     -6/7     -8/7     -10/7    -5/7]
              [                                                ]
              [-1/7    -2/7     -3/7     -4/7     -5/7     -6/7]

> u:=matrix(6,1,[1/49,1/49,1/49,1/49,1/49,1/49]); # THIS IS THE RIGHT
> HAND SIDE VECTOR IN CASE OF F(X)=1

                                  [1/49]
                                  [    ]
                                  [1/49]
                                  [    ]
                                  [1/49]
                             u := [    ]
                                  [1/49]
                                  [    ]
                                  [1/49]
                                  [    ]
                                  [1/49]

> evalm(AINV&*u); # THIS IS THE APPROXIMATE SOLUTION OF THE ORIGINAL
> BOUNDARY VALUE PROBLEM IF F(X)=1

                                [-3 ]
                                [-- ]
                                [49 ]
                                [   ]
                                [-5 ]
                                [-- ]
                                [49 ]
                                [   ]
                                [-6 ]
                                [-- ]
                                [49 ]
                                [   ]
                                [-6 ]
                                [-- ]
                                [49 ]
                                [   ]
                                [-5 ]
                                [-- ]
                                [49 ]
                                [   ]
                                [-3 ]
                                [-- ]
                                [49 ]

> v:=matrix(6,1,[1/343,2/343,3/343,4/343,5/343,6/343]); # THIS IS THE
> RIGHT HAND SIDE VECTOR IN CASE OF F(X)=X

                                  [1/343]
                                  [     ]
                                  [2/343]
                                  [     ]
                                  [3/343]
                             v := [     ]
                                  [4/343]
                                  [     ]
                                  [5/343]
                                  [     ]
                                  [6/343]

> evalm(AINV&*v); # THIS IS THE APPROXIMATE SOLUTION OF THE BOUNDARY
> VALUE PROBLEM IN CASE OF F(X)=X

                                [-8 ]
                                [---]
                                [343]
                                [   ]
                                [-15]
                                [---]
                                [343]
                                [   ]
                                [-20]
                                [---]
                                [343]
                                [   ]
                                [-22]
                                [---]
                                [343]
                                [   ]
                                [-20]
                                [---]
                                [343]
                                [   ]
                                [-13]
                                [---]
                                [343]

> DU:=LUdecomp(A); # THIS IS THE DU FROM THE LDU DECOMPOSITION OF A

                [-2     1       0       0       0       0  ]
                [                                          ]
                [ 0    -3/2     1       0       0       0  ]
                [                                          ]
                [ 0     0      -4/3     1       0       0  ]
          DU := [                                          ]
                [ 0     0       0      -5/4     1       0  ]
                [                                          ]
                [ 0     0       0       0      -6/5     1  ]
                [                                          ]
                [ 0     0       0       0       0      -7/6]

> L:=evalm(A&*inverse(DU)); # THIS IS THE L FROM THE LDU DECOMPOSITION
> OF A.

                [ 1       0       0       0       0      0]
                [                                         ]
                [-1/2     1       0       0       0      0]
                [                                         ]
                [ 0      -2/3     1       0       0      0]
           L := [                                         ]
                [ 0       0      -3/4     1       0      0]
                [                                         ]
                [ 0       0       0      -4/5     1      0]
                [                                         ]
                [ 0       0       0       0      -5/6    1]

> U:=transpose(L); # THIS IS THE U FROM THE LDU DECOMPOSITION OF A. NOTE
> THAT WE HAVE USED PROPOSITION 1.5 SINCE THE A IS SYMMETRIC

                [1    -1/2     0       0       0       0  ]
                [                                         ]
                [0     1      -2/3     0       0       0  ]
                [                                         ]
                [0     0       1      -3/4     0       0  ]
           U := [                                         ]
                [0     0       0       1      -4/5     0  ]
                [                                         ]
                [0     0       0       0       1      -5/6]
                [                                         ]
                [0     0       0       0       0       1  ]

> DI:=evalm(inverse(L)&*A&*inverse(U)); # THIS IS D FROM THE LDU
> DECOMPOSITION OF A

                [-2     0       0       0       0       0  ]
                [                                          ]
                [ 0    -3/2     0       0       0       0  ]
                [                                          ]
                [ 0     0      -4/3     0       0       0  ]
          DI := [                                          ]
                [ 0     0       0      -5/4     0       0  ]
                [                                          ]
                [ 0     0       0       0      -6/5     0  ]
                [                                          ]
                [ 0     0       0       0       0      -7/6]

>