pro gauleg, x1, x2, x, w, n

eps = 3.d-11

x=dblarr(n)
w=dblarr(n)

; Force these to be double precision
  z1=0D & z=0D & xm=0D & xl=0D & pp=0D & p3=0D & p2=0D & p1=0D

  m = (n+1)/2
  xm=0.5D * (x2+x1)
  xl=0.5D * (x2-x1)
  for i=0,m-1 do begin
     z = cos(3.141592654D*(i+0.75D)/(n+0.5D))
     repeat begin
       p1=1D
       p2=0D
       for j=0,n-1 do begin
           p3=p2
           p2=p1
           p1=((2D*j+1D)*z*p2-j*p3)/(j+1)
       end
       pp=n*(z*p1-p2)/(z*z-1D)
       z1=z
       z=z1-p1/pp
     endrep until abs(z-z1) le eps
     x(i)=xm-xl*z
     x(n-1-i)=xm+xl*z
     w(i) = 2D*xl/((1D -z*z)*pp*pp)
     w(n-1-i)=w(i) 
  end
  
end