function xq = gauss_points(Nq,h,xjm1);

  xj = xjm1+h;
  if Nq==2
    xq(1,1) = xjm1 + (h/2)*(1-1/sqrt(3));
    xq(2,1) = xjm1 + (h/2)*(1+1/sqrt(3));
    xq(3,1) =   xj + (h/2)*(1-1/sqrt(3));
    xq(4,1) =   xj + (h/2)*(1+1/sqrt(3));
  else
    xq(1,1) = xjm1 + h*(1/2 - sqrt(15+2*sqrt(30))/(2*sqrt(35)));
    xq(2,1) = xjm1 + h*(1/2 - sqrt(15-2*sqrt(30))/(2*sqrt(35)));
    xq(3,1) = xjm1 + h*(1/2 + sqrt(15-2*sqrt(30))/(2*sqrt(35)));
    xq(4,1) = xjm1 + h*(1/2 + sqrt(15+2*sqrt(30))/(2*sqrt(35)));
    xq(5,1) =   xj + h*(1/2 - sqrt(15+2*sqrt(30))/(2*sqrt(35)));
    xq(6,1) =   xj + h*(1/2 - sqrt(15-2*sqrt(30))/(2*sqrt(35)));
    xq(7,1) =   xj + h*(1/2 + sqrt(15-2*sqrt(30))/(2*sqrt(35)));
    xq(8,1) =   xj + h*(1/2 + sqrt(15+2*sqrt(30))/(2*sqrt(35)));
  end