% QuadgkDemo.m
% Compares |quad| and the newer routine |quadgk| for 
% numerical integration on the interval $[a, b]$. 

clc
% We evaluate
%
% $$ \int_0^1 J_0(x) x e^{x^2} dx$$
%
a = 0;
b = 1;
disp('Evaluate the complicated Bessel integral using quad and quadgk with default tolerance:')
f = @(x) besselj(0,x).*x.*exp(x.^2);
q = quad(f,a,b);
qgk = quadgk(f,a,b);
disp(sprintf('Int(J0(x)*x*exp(x^2),x=%f..%f):',a,b)) 
disp(sprintf('Using quad:   %15.12f\nUsing quadgk: %15.12f',q,qgk)) 
disp('Exact answer:  0.739631779338781') 
disp(' ')

pause
% The |humps| function
%
% $$ \int_{-1}^2 \left( \frac{1}{(x-0.3)^2 + 0.01} + \frac{1}{(x-0.9)^2 + 0.04} - 6 \right) dx $$
%
disp('f(x) = humps(x) on [-1,2]')
f = @(x) humps(x);
a = -1; b = 2;
q = quad(f,a,b);
qgk = quadgk(f,a,b);
disp(sprintf('Using quad:   %15.12f\nUsing quadgk: %15.12f',q,qgk)) 
disp('Exact answer: 26.34496047137833') 
disp(' ')

pause
%
% $$ \int_{0}^{9\pi/2} \cos(x) dx $$
%
disp('f(x) = cos(x) on [0,9*pi/2]')
f = @(x) cos(x);
a = 0; b = 9*pi/2;
tol = 1.e-6;
q = quad(f,a,b,tol);
qgk = quadgk(f,a,b);
disp(sprintf('Using quad:   %15.12f\nUsing quadgk: %15.12f',q,qgk)) 
disp('Exact answer:  1.0') 
disp(' ')

pause
%
% $$ \int_0^1 \sqrt{x} dx $$
%
disp('f(x) = sqrt(x) on [0,1]')
f = @(x) sqrt(x);
a = 0; b = 1;
tol = 1.e-8;
q = quad(f,a,b,tol);
qgk = quadgk(f,a,b);
disp(sprintf('Using quad:   %15.12f\nUsing quadgk: %15.12f',q,qgk)) 
disp('Exact answer:  2/3') 
disp(' ')

pause
%
% $$ \int_0^\pi \left(\tan(\sin(x))-\sin(\tan(x))\right) dx $$
%
disp('f(x) = tan(sin(x))-sin(tan(x)) on [0,pi]')
f = @(x) tan(sin(x))-sin(tan(x));
a = 0; b = pi;
q = quad(f,a,b);
[qgk,errgk] = quadgk(f,a,b);
disp(sprintf('Using quad:   %15.12f\nUsing quadgk: %15.12f',q,qgk)) 
disp('Exact answer: Both Maple and Mathematica give up...') 
disp(sprintf('Error estimate from quadgk: %13.3e',errgk)) 

