function [integrPoints,integrWeights]=computeQuadrature(integrationRule,numIntegrPoints)

%Compute integration points and weights for numerical integration between -1 and 1
if strcmp(integrationRule, 'Gauss-Lobatto')
    if numIntegrPoints==3
        integrPoints=[-1,0,1];
        integrWeights=[1/3,4/3,1/3];
    elseif numIntegrPoints==4
        integrPoints=[-1,-1/sqrt(5),1/sqrt(5),1];
        integrWeights=[1/6,5/6,5/6,1/6];
    elseif numIntegrPoints==5
        integrPoints=[-1,-sqrt(3/7),0,sqrt(3/7),1];
        integrWeights=[1/10,49/90,32/45,49/90,1/10];
    elseif numIntegrPoints==6
        integrPoints=[-1, -0.765055,-0.285232,0.285232, 0.765055, 1];
        integrWeights=[0.066667, 0.378475, 0.554858, 0.554858, 0.378475, 0.066667];
    end
    
% Compute integration points and weights for numerical integration between 0 and 1
    integrPoints=(integrPoints+1)/2;
    integrWeights=integrWeights/2;
    
end