function [integral] = integ_naive (integrand, a, b, N)
% integ_naive calculates an approximation to the integral of
%  the function "integrand" from a to b.
%
%   Input arguments:
%     integrand   name of the function to be integrated in a
%                     string. E.g., integrand='my_function'.  The
%                     function is called by my_function(x).
%     a, b        lower and upper bounds of the integral
%     N           number of subintervals used 
%
%   Output arguments:
%     integral    estimate of the integral

%
% Programmer: Dick Furnstahl  furnstahl.1@osu.edu
%
% Revision history:
%   10-Apr-2005 -- original version for Physics 263
%
% Notes:
%   * The definition of x_pts may be subject to round-off errors
%

% Find Delta x give the number of subintervals N
Delta_x = (b-a) / N;    

% Set up the x points at which the function will be evaluated
%  There will be N+1 points (given N subintervals).
x_pts = a: Delta_x : b;

% Evaluate the function "integrand" at the x points using feval
f_values = feval(integrand, x_pts);   % so f_values is a vector

% Set up the weights for a naive integration 
weights = Delta_x .* ones(1,N+1);  %  start with a vector with N+1 Delta x's 
                                   %  (so this is a 1 x N+1 matrix)
weights(N+1) = 0;       %  set the last weight to zero

% Calculate the estimate to the integral by summing the product of
%  the function and the weights at each of the x points.
integral = sum( f_values .* weights );    % note the .* here