close all; clc; clear; 
randn('seed',0); rand('seed',0);

tol = 1.e-8;

% Problem EPage 45
% 
for n=[3,4,5,6]
  a = rand();
  b = rand();
  t = a + 2*b;
  a = a/(2*t); % scale a and b to make them smaller in magnitude 
  b = b/(2*t); 
  A = diag( a * ones(1,n) ) + diag( b * ones(1,n-1), +1 ) + diag( b * ones(1,n-1), -1 );
  assert( all( abs( eig(A)  ) < 1 ), 'all eigenvalues not less than one' );
  
  % The true value of (I-A)^{-1} 
  % 
  AInv_truth = inv(eye(n) - A);

  % Compute an approximation to (I-A)^{-1} using the power series representation: 
  % 
  k = 1;
  AInv_approx = eye(n);
  delta_AInv = +inf; 
  while delta_AInv > tol
    D_AInv = A^k; % note this is the matrix power
    AInv_approx = AInv_approx + D_AInv; 
    delta_AInv = norm(D_AInv);
    k = k+1;
  end

  disp( sprintf('n= %3d, ||A^(-1) - A^(-1) (Approx)||= %10.6e; with k= %3d terms', n, norm(AInv_truth-AInv_approx), k-1) );
end