function x = QRSolve( A, b)
% solves the system Ax=b by using QR decomposition

[m,n]=size(A);

% auxiliary function: apply Q
function x = ApplyQ( k, alpha, v, x )
    s = 0;
    n = size(v);
    for l = k:n
        s = s + v(l)*x(l);
    end
    s = s / (alpha*v(k));
    for l=k:n
        x(l) = x(l) + s*v(l);
    end
end


% QR decomposition
[A,alpha]=QRDecomposition(A);

% we have to solve the system Rx=Q^Tb, thus we have to compute Q^Tb
for i = 1:n-1
    b = ApplyQ(i, alpha(i), A(:,i), b);
end

[alphaM,alphaN] = size(alpha);

% now we solve the system Rx=Q^Tb
x(m) = b(m) / A(m,m);
for k = m-1:-1:1
    s = b(k);
    for j = k+1:n
        s = s - A(k,j)*x(j);
    end
    x(k) = s/alpha(k);
end

end

