h = 0.001;
L = 100;                            % exp(-100) is really small
x = (0:h:L)';

ex = exp(-x);
A  = [ex, x.*ex, x.^2.*ex, x.^3.*ex];

% modify for trapezoid rule estimate of inner product
At      = A; 
At(1,:) = sqrt(.5)*At(1,:);
At(end,:) = sqrt(.5)*At(end,:);
At      = sqrt(h)*At;

% check value of some column norms
% compare [exact  trapezoid  left-hand] for ||a1||^2 and ||a2||^2:
% [0.5,  At(:,1)'*At(:,1),  A(:,1)'*A(:,1)*h]
% [0.25, At(:,2)'*At(:,2),  A(:,2)'*A(:,2)*h]

[Q,R] = qr(At,0);                   % *reduced* QR, naturally

% change signs to match Odile's answer, which has r_ij >= 0:
S     = diag(sign(diag(R)));        % S = S^-1 and S is unitary
Q     = Q*S;
R     = S*R;

disp("Odile's by-hand answer is 1/sqrt(2) times this matrix:")
format rat                          % show with fractions
Rfract = [1 0.5 0.5 0.75; 0 0.5 1 2.25; 0 0 0.5 2.25; 0 0 0 .75]
format                              % restore default format
Rexact = (1/sqrt(2)) * Rfract;

disp(["distance between numerical R (trapezoid and built-in QR) and by-hand R: " ...
     num2str(norm(R-Rexact))])

if norm(R-Rexact) < 1.0e-4
  disp("both numerical and by-hand are correct with very high probability!")
end