function Q = dblquad3(intfcn,xmin,xmax,ymin,ymax,zmin,zmax,tol,quadf1,quadf2,varargin)
if nargin < 7, error(’Requires at least five inputs’); end
if nargin < 8 | isempty(tol), tol = 1.e-6; end
if nargin < 9 | isempty(quadf1), quadf1=@quad; end
if nargin < 10 | isempty(quadf2), quadf2=@dblquad; end
intfcn = fcnchk(intfcn);
trace = [];
Q = feval(quadf1,@innerintegral1,zmin,zmax,tol,trace,intfcn, ...
xmin,xmax,ymin,ymax,tol,quadf2,varargin{:);
%---------------------------------------------------------------------------
function Q = innerintegral1(z,intfcn,xmin,xmax,ymin,ymax,tol,quadf2,varargin)
Q = zeros(size(z));
for i = 1:length(z)
Q(i) = feval(quadf2,intfcn, xmin,xmax,ymin,ymax,tol,@quad,z(i),varargin{:);
end
function Q = dblquad4(intfcn,xmin,xmax,ymin,ymax,zmin,zmax,rmin,rmax,tol,quadf3,quadf4,varargin)
if nargin < 9, error(’Requires at least five inputs’); end
if nargin < 10 | isempty(tol), tol = 1.e-6; end
if nargin < 11 | isempty(quadf3),quadf3=@quad; end
if nargin < 12 | isempty(quadf4),quadf4=@dblquad3; end
intfcn = fcnchk(intfcn);
trace = [];
Q = feval(quadf3,@innerintegral2,rmin,rmax,tol,trace,intfcn, ...
xmin,xmax,ymin,ymax,zmin,zmax,tol,quadf4,varargin{:);
%---------------------------------------------------------------------------
function Q = innerintegral2(r,intfcn,xmin,xmax,ymin,ymax,zmin,zmax,tol,quadf4,varargin)
Q = zeros(size(r));
for i = 1:length(r);
Q(i) = feval(quadf4,intfcn,xmin,xmax,ymin,ymax,zmin,zmax,tol,@quad,@dblquad,r(i),varargin{:);
end;
