   function [xi,wi]=gauleg(a,b,n)
%  Gauss-Legndre quadrature, w(x)=1 on range  (-1,1).
   m=(n+1)/2 ;
   xm=(a+b)/2;
   xl=(b-a)/2;
   EPS = 3e-14;
   for i=1:m
       z=cos(pi*(i-.25)/(n+0.5));
       z1 = z+10;
       while (abs(z-z1)> EPS) 
           p1=1 ;
           p2=0 ;
           for j=1:n
               p3=p2;
               p2=p1;
               p1=((2*j-1)*z*p2-(j-1)*p3)/j ;
           end
           pp=n*(z*p1-p2)/(z^2-1);
           z1=z;
           z=z1-p1/pp;
       end
       xi(i)=xm-xl*z ;
       xi(n+1-i)=xm+xl*z ;
       wi(i)=2.d0*xl/((1.d0-z*z)*pp*pp);
       wi(n+1-i)=wi(i) ;
   end 
