% testIntegrateFunc.m
% Tests integrateFunc by comparing its output with that of integral
% J.M.Soler. Nov.2011

% Set integration interval between zero and first root of j1(x)
clear all
x1 = 0;
x2 = bisectionRoot( @j1, pi, 2*pi, 1.e-12 );

% Find integral using integrateFunc
dx = (x2-x1)/100;
za = integrateFunc( @j1, x1, x2, dx )

% Find integral using function integral
x = x1:dx:x2;           % mesh of integration points
y = j1(x);              % values of j1 at the mesh points
z = integral( x, y );   % values of the integral at mesh points
zb = z(end)             % value of z at last point
