% Use 3pt Gauss quadrature to integrate function f(x) to 5th order
% This way we compute finite-volume representations of a given function f(x)
function w = FVQuad(f,x,h) % f(x) must be vectorized!
  
  if(1) % Use a finite-volume interpretation  
    w = GaussPoint(0,8/9,f,x,h);
    w = w + GaussPoint(+sqrt(3/5),5/9,f,x,h);
    w = w + GaussPoint(-sqrt(3/5),5/9,f,x,h);
    w = w/h; % The value is the average not the integral
  else
    w = f(x); % Simple finite-difference pointwise interpretation
  end  
    
end

function dw = GaussPoint(gauss,weight,f,x,h)
    points=((x-h/2)*(1-gauss) + (x+h/2)*(1+gauss))/2;
    dw=weight*f(points)*h/2;
end
