n = 10; 
    A = DiagonalMatrix[Table[2, {i, 1, n - 1}]];
For[i = 1, 
        i < n - 1, i++, 
A[[i, i + 1]] = -1; 
        A[[i + 1, 
              i]] = -1;
];
B = Inverse[A];
Print[MatrixForm[A]];
Print["Inverse of Sparse is Full"];
Print[MatrixForm[B]];

h = n^2;
Print["Eigenvalues of A ~ Eigenvalues of -D2"];
Sort[evals = Eigenvalues[h A] // N]
Table[i^2  Pi^2, {i, 1, n - 1}] // N

ev[[8]]

Print["Eigenvectors of A ~ Eigenfunctions of -D2"];
evects = Eigenvectors[h A] // N;
For[i = 1, 
        i ≤ n - 1, i++, 
norm = 
          Sqrt[evects[[i]] . 
              evects[[i]]]; 
evects[[i]] /= 
          norm;
];
e8 = Join[{0}, evects[[8]], {0}];
e8x = 
        Table[{i/n, e8[[i + 1]]}, {i, 0, 
            n}];
Print[{n, Length[e8x]}];
ListPlot[e8x, PlotJoined -> True];

Print["Symmetric Matrix has orthonormal eigenvectors"];
Print[MatrixForm[Chop[evects . Transpose[evects]]]];

Print["Solving system approximates 2nd order BVP solution"];
b = Table[1, {i, 1, n - 1}];
c = Join[{0}, x = B . b, {0}];
cc = 
        Table[{i/n, c[[i + 1]]}, {i, 0, n}];
ListPlot[cc, PlotJoined -> True];

Print["More efficient to solve Sparse system directly"];
m = n - 1;
DD = Table[2, {i, 1, m}];
UU = Table[-1, {i, 1, m}];
LL = UU;
bb = b;
For[i = 2, 
        i ≤ 
          m, i++, 
fact = -LL[[i]]/
            DD[[i - 1]]; 
DD[[i]] += 
          fact*UU[[i - 
                1]]; 
bb[[i]] += 
          fact*bb[[i - 
                1]]; 
];

\(\(bb[\([m]\)] /= DD[\([m]\)];\)
\)
For[i = m - 1, 
        i ≥ 
          1, i--, 
bb[[i]] = (bb[[i]] - UU[[i]]*bb[[i + 1]])/
            DD[[i]]; 
];
Print["Product A^-1 b:"];
Print[x];
Print["Solving system A.x = b"];
Print[bb];

f[{x_, y_}] := Sin[x] Sin[y];
Print["Find min/max/saddle points of f"];
Plot3D[f[{x, y}], {x, -Pi, Pi}, {y, -Pi, Pi}]

Clear[g, h];
g[{x_, y_}] := {D[f[{s, t}], s], D[f[{s, t}], t]}  /.  
            s -> x  /.  t -> y;
h[{x_, y_}] := {D[g[{s, t}], s], D[g[{s, t}], t]}  /.  
            s -> x  /.  t -> y;

{x, y} = {0.1, 0.1};

{x, y} = {1, 1};

{x, y} = {-1, 1};

go = 1; tol = 10^-10; cnt = 0; 
    maxits = 10;
gg = g[{x, y}] // N;
hh = h[{x, y}] // N;
\(\(While[
        go == 
          1, 
\(cnt++\); 

 (*\ 
          BETTER\ to\ solve\ system\ *) 
{x, y} -= 
            Inverse[hh] . gg // N; 

gg = 
          g[{x, y}] // N; 
hh = 
          h[{x, y}] // N; 

norm\  = \ 
          Sqrt[gg . gg]; 
ev\  = \ 
          Eigenvalues[
            hh]; 

Print["cnt = ", 
          cnt, ", norm = ", \ norm, \ ", ev = ", 
          ev]; 
If[\((cnt ≥ 
                maxits)\) || \((norm < tol)\), go = 0, 
          Null];

];\)\)