function [nodes,weights] = GaussLaguerreFormula(n_nodes)

%

Pols = LaguerrePolynomials(n_nodes + 1);
Pol = Pols{end - 1};
nodes = sort(roots(Pol))';
Pol_weights = Pols{end};
weights = nodes ./ (((n_nodes + 1) ^ 2) * ((polyval(Pol_weights,nodes)) .^ 2));

end