function [r, w] = Gauss(l)
% RETURN COORDINATES OF GAUSS QUADRATURE POINTS AND THEIR WEIGHTS
%
% SYNTAX: [r, w] = Gauss(l)
%
% INPUT PARAMETERS
%    l - Order of Gauss quadrature ( 1 =< r <= 10 )
%
% OUTPUT PARAMETERS
%    r - Coordinates of Gauss quadrature points
%    w - Weights at Gauss quadrature points
%
xg = [
               0.0              0.0              0.0              0.0              0.0              0.0              0.0              0.0              0.0              0.0
  -0.5773502691896  0.5773502691896              0.0              0.0              0.0	            0.0              0.0              0.0              0.0              0.0
  -0.7745966692415  0.0000000000000  0.7745966692415              0.0              0.0	            0.0              0.0              0.0              0.0              0.0 
  -0.8611363115941 -0.3399810435849  0.3399810435849  0.8611363115941              0.0	            0.0              0.0              0.0              0.0              0.0
  -0.9061798459387 -0.5384693101057              0.0  0.5384693101057  0.9061798459387              0.0              0.0              0.0              0.0              0.0 
  -0.9324695142032 -0.6612093864663 -0.2386191860832  0.2386191860832  0.6612093864663  0.9324695142032              0.0              0.0              0.0              0.0 
  -0.9491079123428 -0.7415311855994 -0.4058451513774              0.0  0.4058451514774	0.7415311855994  0.9491079123428              0.0              0.0              0.0
  -0.9602898564975 -0.7966664774136 -0.5255324099163 -0.1834346424957  0.1834346424957	0.5255324099163  0.7966664774136  0.9602898564975              0.0              0.0 
  -0.9681602395076 -0.8360311073266 -0.6133714327006 -0.3242534234038              0.0	0.3242534234038  0.6133714327006  0.8360311073266  0.9681602395076              0.0
  -0.9739065285172 -0.8650633666890 -0.6794095682990 -0.4333953941292 -0.1488743389816  0.1488743389816  0.4333953941292  0.6794095682990  0.8650633666890  0.9739065285172   
      ];

wg  = [
              2.0,             0.0,             0.0,             0.0,              0.0,             0.0,             0.0,             0.0, 			  0.0,             0.0  
              1.0,             1.0,             0.0,             0.0,              0.0,             0.0,             0.0,             0.0, 			  0.0,             0.0  
  0.5555555555556, 0.8888888888889, 0.5555555555556,             0.0,              0.0,             0.0,             0.0,             0.0,			  0.0,             0.0  
  0.3478548451375, 0.6521451548625, 0.6521451548625, 0.3478548451375,              0.0,             0.0,             0.0,             0.0,			  0.0,             0.0  
  0.2369268850562, 0.4786286704994, 0.5688888888889, 0.4786286704994,  0.2369268850562,             0.0,             0.0,             0.0,			  0.0,             0.0  
  0.1713244923792, 0.3607615730481, 0.4679139345727, 0.4679139345727,  0.3607615730481, 0.1713244923792,             0.0,             0.0,			  0.0,             0.0  
  0.1294849661689, 0.2797053914893, 0.3818300505051, 0.4179591836735,  0.3818300505051, 0.2797053914893, 0.1294849661689,             0.0,			  0.0,             0.0  
  0.1012285362904, 0.2223810344534, 0.3137066458779, 0.3626837833784,  0.3626837833784, 0.3137066458779, 0.2223810344534, 0.1012285362904,			  0.0,             0.0  
  0.0812743883616, 0.1806481606949, 0.2606106964029, 0.3123470770400,  0.3302393550013, 0.3123470770400, 0.2606106964029, 0.1806481606949, 0.0812743883616,             0.0  
  0.0666713443087, 0.1494513491506, 0.2190863625160, 0.2692667193100,  0.2955242247148, 0.2955242247148, 0.2692667193100, 0.2190863625160, 0.1494513491506, 0.0666713443087  
      ];

if (l > 10)
   l = 10;
end
   
if (l < 1)
   l = 1;
end

r = zeros(1,l);
r(1,1:l) = xg(l,1:l);
w(1,1:l) = wg(l,1:l);