function [W,Q]=quadrature(qord,dim)
qpoints=zeros(qord^dim,dim);
qweights=zeros(qord^dim,1);
r1pt=zeros(qord,1);
r1wt=zeros(qord,1);

switch (qord)
    case 1
        r1pt(1)=0.0000000000;
        r1wt(1)=2.0000000000;
    case 2
        r1pt(1)=0.57735;
        r1pt(2)=-0.57735;
        
        r1wt(1)=1.00000;
        r1wt(2)=1.00000;
    case 3
        r1pt(1)=0.77459;
        r1pt(2)=-0.77459;
        r1pt(3)=0.00000;
        
        r1wt(1)=0.55555;
        r1wt(2)=0.55555;
        r1wt(3)=0.88888;
end
n=1;
if(dim==1)
    for i=1:qord
    qpoints(n,:)=[r1pt(i)];
    qweights(n)=r1wt(i);
    n=n+1;
    end
elseif(dim==2)
    for i=1:qord
        for j=1:qord
            qpoints(n,:)=[r1pt(i),r1pt(j)];
            qweights(n)=r1wt(i)*r1wt(i);
            n=n+1;
        end
    end
end
Q=qpoints;
W=qweigths;
        