function [nodes,weights] = GaussLobattoFormula(n_nodes)

% computes nodes and corresponding weights for the Gauss-Lobatto guadrature, when number of nodes is given = n_nodes

deg = n_nodes - 1;
Pol = LegendrePolynomial(deg);
dPol = polyder(Pol);
roots_dPol = sort(roots(dPol))';
nodes = [-1 roots_dPol 1];
weights = 2 ./ (deg * (deg + 1) .* polyval(Pol,nodes) .^ 2);

end
