function a = myqrfit(x, y)
% MYQRFIT - least squares fitof noisy (x,y) data using QR factorization
%          Return the coefficients a of the expansion 
%          f(x) = a(1)./x + a(2)*cos(x) + a(3)*exp(-x.*x/4)
   A = [1 ./(x') cos(x') exp(-x'.*x'/4)];
   [Q, R] = myqrfact(A);
   d = Q'*y';
   a = R\d;
end
