clear all
n    = 30;
x    = zeros(n,1);
x(1) = 1;
x(2) = 1/3;
for ii = 3:n
    x(ii)   = 13 / 3 * x(ii-1) - 4 / 3 * x(ii-2);
    x_ast   = (1/3)^(ii-1);
    Rel_err_x = abs(x_ast-x(ii)) / x_ast;
  fprintf('x(%2.0f) = %16.8d, x*(%2.0f) = %16.8d, RelErr(%2.0f) = %10.4d \n'...
          ,ii,x(ii),ii,x_ast,ii,Rel_err_x);
end