From CFD-Wiki
function [ xa,ya,wt,nq ] = GQuad2( ord )
%GQuad2
% Returns arrays of integration points and weights for specified degree
% Usage:
% [xa,ya,wt,nq] = GQuad2( ord )
% where
% ord - highest degree of term integrated exactly on interval [-1,1]
%
if ord<1, error('integration order < 1');
elseif ord>15, error('max quad data order = 15');
end
nt = fix(ord/2)+1;
switch nt;
case {1, 2} % degree 3
Aq=[-0.5773502691896257, 0.5773502691896257];
Hq=[ 1.0, 1.0];
case 3 % degree 5
Aq=[-0.7745966692414833, 0.000000000000000, 0.7745966692414833];
Hq=[ 0.5555555555555558, 0.8888888888888888,0.5555555555555558];
case 4 % degree 7
Aq=[-0.8611363115940526,-0.3399810435848563, 0.3399810435848563, ...
0.8611363115940526];
Hq=[ 0.3478548451374538, 0.6521451548625461, 0.6521451548625461, ...
0.3478548451374538];
case 5 % degree 9
Aq=[-0.9061798459386640,-0.5384693101056831, 0.0000000000000000, ...
0.5384693101056831, 0.9061798459386640];
Hq=[ .23692688505618909, .47862867049936647, .56888888888888889, ...
.47862867049936647, .23692688505618909];
case 6 % degree 11
Aq=[-.932469514203152,-.661209386466265,-.238619186083197, ...
.238619186083197, .661209386466265, .932469514203152];
Hq=[ .171324492379170, .360761573048139, .467913934572691, ...
.467913934572691, .360761573048139, .171324492379170];
case 7 % degree 13
Aq=[-0.9491079123427585,-0.7415311855993945,-0.4058451513773972, ...
0.0000000000000000, 0.4058451513773972, 0.7415311855993945, ...
0.9491079123427585];
Hq=[ 0.1294849661688699, 0.2797053914892767, 0.3818300505051190, ...
0.4179591836734694, 0.3818300505051190, 0.2797053914892767, ...
0.1294849661688699];
case 8 % degree 15
Aq=[-0.9602898564975363,-0.7966664774136267,-0.5255324099163290, ...
-0.1834346424956498, 0.1834346424956498, 0.5255324099163290, ...
0.7966664774136267, 0.9602898564975363];
Hq=[ 0.1012285362903761, 0.2223810344533745, 0.3137066458778873, ...
0.3626837833783620, 0.3626837833783620, 0.3137066458778873, ...
0.2223810344533745, 0.1012285362903761];
end
nt=size(Aq,2);
xa = repmat(Aq,nt,1);
ya = xa';
wt = repmat(Hq,nt,1);
wt = wt.*wt';
nq = nt*nt;
return;
end