matlab quad 函数代码中的y(1) , y(7)什么意思, 跟数值积分分成的n个节点有什么关系，请高手帮忙指点

function [Q,fcnt] = quad(funfcn,a,b,tol,trace,varargin)
%QUAD Numerically evaluate integral, adaptive Simpson quadrature.
% Q = QUAD(FUN,A,B) tries to approximate the integral of scalar-valued
% function FUN from A to B to within an error of 1.e-6 using recursive
% adaptive Simpson quadrature. FUN is a function handle. The function
% Y=FUN(X) should accept a vector argument X and return a vector result
% Y, the integrand evaluated at each element of X.
%
% Q = QUAD(FUN,A,B,TOL) uses an absolute error tolerance of TOL
% instead of the default, which is 1.e-6. Larger values of TOL
% result in fewer function evaluations and faster computation,
% but less accurate results. The QUAD function in MATLAB 5.3 used
% a less reliable algorithm and a default tolerance of 1.e-3.
%
% Q = QUAD(FUN,A,B,TOL,TRACE) with non-zero TRACE shows the values
% of [fcnt a b-a Q] during the recursion. Use [] as a placeholder to
% obtain the default value of TOL.
%
% [Q,FCNT] = QUAD(...) returns the number of function evaluations.
%
% Use array operators .*, ./ and .^ in the definition of FUN
% so that it can be evaluated with a vector argument.
%
% Notes:
% Function QUADL may be more efficient with high accuracies and smooth
% integrands.
% Function QUADV vectorizes QUAD for array-valued FUN.
%
% Example:
% Q = quad(@myfun,0,2);
% where myfun.m is the M-file function:
% %-------------------%
% function y = myfun(x)
% y = 1./(x.^3-2*x-5);
% %-------------------%
%
% or, use a parameter for the constant:
% Q = quad(@(x)myfun2(x,5),0,2);
% where myfun2 is the M-file function:
% %----------------------%
% function y = myfun2(x,c)
% y = 1./(x.^3-2*x-c);
% %----------------------%
%
% Class support for inputs A, B, and the output of FUN:
% float: double, single
%
% See also QUADV, QUADL, DBLQUAD, TRIPLEQUAD, TRAPZ, FUNCTION_HANDLE.

% Based on "adaptsim" by Walter Gander.
% http://www.inf.ethz.ch/personal/gander
%
% Reference:
% [1] W. Gander and W. Gautschi, Adaptive Quadrature - Revisited,
% BIT Vol. 40, No. 1, March 2000, pp. 84-101.
%
% Copyright 1984-2006 The MathWorks, Inc.
% $Revision: 5.26.4.7 $ $Date: 2006/04/03 17:10:41 $

f = fcnchk(funfcn);
if nargin < 4 || isempty(tol), tol = 1.e-6; end;
if nargin < 5 || isempty(trace), trace = 0; end;
if ~isscalar(a) || ~isscalar(b)
error('MATLAB:quad:scalarLimits',...
'The limits of integration must be scalars.');
end

% Initialize with three unequal subintervals.
h = 0.13579*(b-a);
x = [a a+h a+2*h (a+b)/2 b-2*h b-h b];
y = f(x, varargin{:});
fcnt = 7;

% Fudge endpoints to avoid infinities.
if ~isfinite(y(1))
y(1) = f(a+eps(superiorfloat(a,b))*(b-a),varargin{:});
fcnt = fcnt+1;
end
if ~isfinite(y(7))
y(7) = f(b-eps(superiorfloat(a,b))*(b-a),varargin{:});
fcnt = fcnt+1;
end

这个是按着辛普森公式在计算，辛普森公式跟节点有关的，Y1应该是Y的第一个值，就是说第一次算出来的值，Y7是Y积分结果矩阵的第七个值，这个值都跟节点有关系。不好意思。这个辛普森公式在计算方法（或者数值分析里有详细介绍），我给忘记了！所以只能回答这么多！

应这样改： 一个文件中内容如下： function y = fun(x) y = x.^2+3*x.^3-1; end 另一个主文件内容如下： quad('fun',-1,1)