
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 = abs(x_ast-x(ii)) / x_ast;
    fprintf('x(%2.0f) = %20.8d,   x^\ast(%2.0f) = %20.8d,  RelErr(%2.0f) = %14.4d \n', ...
        ii,x(ii),ii,x_ast,ii,Rel_err);
end
