quadshow(fun, a, b, varargin)

function Q = quadshow(fun, a, b, varargin)
% QUADSHOW Display the points used to integrate a function using QUAD
%    Q = QUADSHOW(FUN, A, B) uses QUAD to integrate the function FUN from A
%    to B. While the integration is being performed, QUADSHOW displays a
%    graph of the function indicating the points at which QUAD is
%    integrating the function.  The function FUN must accept a vector
%    argument X and return a vector result Y that is the same size as X,
%    exactly as required by QUAD.
%
%    Q = QUADSHOW(FUN, A, B, TOL) is the same as Q = QUAD(FUN, A, B, TOL)
%    except that QUADSHOW displays the function as described above.
%
%    Q = QUADSHOW(FUN, A, B, TOL, TRACE) is the same as
%    Q = QUAD(FUN, A, B, TOL, TRACE) except that QUADSHOW displays the
%    function as described above.
%
%    Examples:
%
%    % Example 1
%    Q = quadshow(@(x) x.^2, 0, 1);
%
%    % Example 2
%    figure;
%    subplot(3, 1, 1);
%    Q1 = quadshow(@(x) (x.^2).*exp(-x.^2), -4, 4)
%    subplot(3, 1, 2);
%    Q2 = quadshow(@(x) (x.^2).*exp(-x.^2), -10, 10)
%    subplot(3, 1, 3);
%    Q3 = quadshow(@(x) (x.^2).*exp(-x.^2), -1000, 1000)
%
% See Also: QUAD
%
%    Note: I could have modified this to allow users to pass an input
%    argument to select which quadrature routine (QUAD or QUADL) to use to
%    integrate the function, but I chose to keep the syntax for QUADSHOW as
%    close to the syntax of QUAD as I could.

% Copyright 2007 The MathWorks, Inc.

% Create the line object that will display the function
h = line(NaN, NaN, 'Marker', 'o');

% The axes only need to cover the X range between the two limits
set(gca, 'XLim', [a b]);

% Title the axes
t1 = title(sprintf('Integral of %s from %g to %g', func2str(fun), a, b), ...
    'Interpreter', 'none');

% Perform the actual quadrature using a nested function that will call the
% user-specified function
Q = quad(@nestfun, a, b, varargin{:});

% Update title
set(t1, 'String', sprintf('Integral of %s from %g to %g = %g', func2str(fun), a, b, Q));

% Begin nested function nestfun
    function z = nestfun(t)
        
        % Evalute the function.  I'll try to be nice and catch two of the
        % common QUAD problems right here.
        try
            z = fun(t);
        catch
            S = lasterror;
            % Common problem 1:  If you used ^ instead of .^ or * instead
            % of .*, suggest the elementwise operators in the error
            % message.  Compare @(x) x^2 and @(x) x.^2
            if ismember(S.identifier, {'MATLAB:square', 'MATLAB:innerdim'})
                S.message = sprintf('%s\n%s\n\n%s', ...
                    ['The integrand function reported an error.  If your', ...
                    ' function uses ^ or *, it may need to'], ...
                    ['use .^ or .* to perform elementwise operations', ...
                    ' instead.  The error the integrand reported was: '], ...
                    S.message);
            end
            rethrow(S);
        end
        % Common problem 2:  QUAD requires the function you pass to it to
        % return output the same size as the input.  This commonly comes up
        % with something like @(x) 1
        if ~isequal(size(t), size(z))
            error('quadshow:OutputSameSizeAsInput', ...
                ['The output of the function you passed to QUADSHOW ', ...
                'should be the same size as its input.']);
        end
        
        % This TRY/CATCH block is in place in case we try to manipulate the
        % line, but it's already been deleted
        try
            % Append the new X and Y coordinates to the existing line data
            XD = [get(h, 'XData') t];
            YD = [get(h, 'YData') z];
            % Sort the XData so we don't get zig-zags
            [sx, ord] = sort(XD);
            % Sort the YData the same way
            sy = YD(ord);
            set(h, 'XData', sx, 'YData', sy);
            % Pause to allow the user to see the new point; adjust this to
            % speed up or slow down the display
            pause(0.5)
        catch
        end
    end
% End nested function nestfun
end