function orthonormal = PolynomialGramSchmidtOrthonormalization(polynomials)

% vectors = matrix with polynomials (vectors) in its rows to be orthonormalized

n = size(polynomials,2);

orthogonal = zeros(n);
orthonormal = orthogonal;

orthogonal(1,:) = polynomials(1,:);
for i = 2:n
    temp = polynomials(i,:);
    for j = 1:(i - 1)
        temp = temp - PolynomialOrthogonalProjection(polynomials(i,:),orthogonal(j,:));
    end
    orthogonal(i,:) = temp;
end

for i = 1:n
    orthonormal(i,:) = orthogonal(i,:) / PolynomialNorm(orthogonal(i,:),[-1, 1]);
end

end