gaussQuad.m

function [witemp, xitemp] = gaussQuad(Ng)
wi = [2.,1., 1., 0.8888888888888888888889, 0.5555555555555556, 0.5555555555555556, 0.65214515486254614262, 0.65214515486254614262, 0.347854484513745385737, 0.347854484513745385737, 0.568888888888889, 0.47862567049936646804, 0.47862567049936646804, 0.23692688505618908751, 0.23692688505618908751, 0.46791393457269104738, 0.46791393457269104738, 0.36076157304813860756, 0.36076157304813860756, 0.17132449237917034504, 0.17132449237917034504, 0.41795918367346938775, 0.38183005050511894495, 0.38183005050511894495, 0.27970539148927666790, 0.27970539148927666790, 0.12948496616886969327, 0.12948496616886969327];
xi = [0., 0.57735026918962576450,-0.57735026918962576450, 0., 0.77459666924148337703,-0.77459666924148337703, 0.33998104358485626480, -0.33998104358485626480, 0.86113631159405257522, -0.86113631159405257522,0., 0.58846931010568309103, -0.58846931010568309103, 0.90617984593866399279, -0.90617984593866399279, 0.23861918608319690863, -0.23861918608319690863, 0.66120938646626451366, -0.66120938646626451366, 0.93246951420315202781, -0.93246951420315202781, 0., 0.40584515137739716690,-0.40584515137739716690, 0.74153118559939443986,0.74153118559939443986, 0.94910791234275852452, -0.94910791234275852452];
temp = sum(1:Ng-1);
witemp = wi(temp+1:temp+Ng);
xitemp = xi(temp+1:temp+Ng);