%
% this code demonstrates how the different versions of 
% Simpson's rule work on example functions
%
% Written by Ming Gu for Math 128a, Fall 2008
%
format short e;
interv = [0;1];
tol0   = 1e-8;
%
% use matlab function quad to evalute integral
%
Int.q  = quad(@(x)xlnx(x),interv(1),interv(2),tol0);
ff = ['Integrating over interval [',num2str(interv(1)), ', ', num2str(interv(2)), ']'];
disp(ff);
disp([]);
%
% Simpson's rule
%
Int.simpson = Simpson(@(x)xlnx(x),interv);
ff1 = ['Simpson rule =  ',num2str(Int.simpson)];
ff2 = ['Simpson rule error = ', num2str(Int.simpson-Int.q)];
ff  = strvcat(ff1,ff2);
disp(ff);
%
% Composite Simpson's rule
%
n           = 501;
Int.comp    = CompSimpson(@(x)xlnx(x),n,interv);
ff1 = ['Composite Simpson rule =  ',num2str(Int.comp), ' with ', num2str(n), ' points'];
ff2 = ['Composite Simpson rule error = ', num2str(Int.comp-Int.q)];
ff  = strvcat(ff1,ff2);
disp(ff);
%
% Adaptive Simpson's rule
%
tol         = 5e-7;
L           = 20; 
[Int.adap,flg,fcnt,level] =AdaptSimpson(@(x)xlnx(x),intv,tol,L);
ff1 = ['Adaptive Simpson rule =  ',num2str(Int.adap), ' with tol = ', num2str(tol)];
ff2 = ['Adaptive Simpson evaluated function at ', num2str(fcnt), ' points on ', num2str(level) ' levels'];
ff3 = ['Adaptive Simpson error =  ',num2str(Int.adap-Int.q)];
ff  = strvcat(ff1,ff2,ff3);
disp(ff);
%
% Int is the data structure containing all the integral values
%