function Itpc = quadc(a,b,M,f,fp,choice,varargin)
% Choice: 1, trapezpoidal, order 2
% 2, Gauss composite order 4
% 3, "" """ order 6
% 4, Simpson, order 4
H = (b-a)/M;
switch choice
case 1,
% Trapezoidal, interpolatory order 1
% Nodes: 0, 1
% Weights: 1/2, 1/2
x = linspace(a,b,M+1);
fpm = feval(f,x,varargin{:});
fpm(2:end-1) = 2*fpm(2:end-1);
Itpc = 0.5*H*sum(fpm);
error = H^2*fp*(b-a)/12
case 2,
% Gauss-Legendre, order 3
% Nodes: zeros of P2 = x^2-1/3
% Weights: 1/2, 1/2
z = linspace(a,b,M+1)
x = z(1:end-1)+H/2*(1-1/sqrt(3))
x = [x,x+H/sqrt(3)]
fpm = feval(f,x,varargin{:});
Itpc = 0.5*H*sum(fpm);
error = H^4*fp*(b-a)/4320
case 3
% Gauss-Legendre, order 5
% Nodes: zeros of P3 = x^3-x*3/5
% Weights: 5/18, 8/18, 5/18
z = linspace(a,b,M+1)
x = z(1:end-1)+H/2*(1-sqrt(3/5))
x = [x,x+H*sqrt(3/5)]
fpm1 = feval(f,x,varargin{:});
Itpc1 = 5*sum(fpm1);
y = z(1:end-1)+H/2
fpm2 = feval(f,y,varargin{:});
Itpc2 = 8*sum(fpm2);
Itpc = (Itpc1+Itpc2)*H/18;
case 4
% Simpson's, interpolatory order 3 (accurate order 4)
% Nodes: 0, 1/2, 1
% weights: 1/6, 4/6, 1/6
x = linspace(a,b,M+1);
fpm = feval(f,x,varargin{:});
fpm(2:end-1) = 2*fpm(2:end-1);
Isic = H*sum(fpm)/6;
x = linspace(a+H/2,b-H/2,M);
fpm = feval(f,x,varargin{:});
Itpc = Isic + 2*H*sum(fpm)/3;
error = -(1/16)*H^4*fp*(b-a)/180
otherwise
disp('no such case');
end