function [nodes, weights] = legendre_quad(lo, hi, num_nodes)

  
%  usage:  [nodes, weights] = legendre_quad(lo,hi,num_nodes)
%
%  Calculates the Legendre polynomial nodes and weights 
%  for given range and num_nodes.
%  From NR p. 152
  
  tol = (eps);
  m=(num_nodes+1)/2;
  xm=0.5*(hi+lo);
  xl=0.5*(hi-lo);
  
  for i=1:m
  
    z=cos(pi*(i-0.25)/(num_nodes+0.5)); 
     % Initial guess at z
    z1 = z+.1;
    while abs(z-z1) > tol                 
      % then refine Newton's method...
      p1=1.0;
      p2=0.0;
      
      for j=1:num_nodes
     	p3=p2;
	    p2=p1;
	    p1=((2*j-1)*z*p2-(j-1)*p3)/j;
      end
      
      pp=num_nodes*(z*p1-p2)/(z*z-1.0);
      z1=z;
      z=z1-p1/pp;
      
    end
      
    nodes(i)=xm-xl*z;
    nodes(num_nodes+1-i)=xm+xl*z;
    weights(i)=2*xl/((1-z*z)*pp*pp);
    weights(num_nodes+1-i)=weights(i);
 
end
nodes=nodes';
weights=weights';
