% exercise 9.1 (c), Trefethen & Bau, page 64
nulist = 1:15;
dx = zeros(size(nulist));  err = dx;  % space for dx, errors

for j = 1:length(nulist);
  dx(j) = 2^(-nulist(j));   disp(["dx = " num2str(dx(j)) ":"])
  x = (-1.0:dx(j):1.0)';
  A = [x.^0 x.^1 x.^2 x.^3];
  [Q,R] = qr(A,0);
  Papprox = Q(:,4) / Q(end,4);    Pexact = 2.5 * x.^3 - 1.5 * x;
  err(j) = max(abs(Papprox - Pexact));
end

loglog(dx,err,'o',"markersize",2)