%% Chapter 5
%% Piecewise Interpolation and Calculus
% _This is a MATLAB Live Editor version of the *Fundamentals of Numerical Computation* 
% Jupyter Notebooks found at_ <https://tobydriscoll.net/project/fnc/ _https://tobydriscoll.net/project/fnc/_>_._ 
% Polynomials
% 
% Example 5.1.1
% Here are some points that we could consider to be observations of an unknown 
% function on $[-1,1]$.

clear
n = 5;
t = linspace(-1,1,n+1)';  y = t.^2 + t + 0.05*sin(20*t);
plot(t,y,'.')
%% 
% The polynomial interpolant, as computed using |polyfit|, looks very sensible. 
% It's the kind of function you'd take home to meet your parents. The interpolant 
% plot is created using |fplot|.

c = polyfit(t,y,n);     % polynomial coefficients
p = @(x) polyval(c,x);
hold on, fplot(p,[-1 1])
%% 
% But now consider a different set of points generated in almost exactly the 
% same way.

n = 18;
t = linspace(-1,1,n+1)';  y = t.^2 + t + 0.05*sin(20*t);
clf, plot(t,y,'.')
%% 
% The points themselves are unremarkable. But take a look at what happens to 
% the polynomial interpolant.

c = polyfit(t,y,n);     % polynomial coefficients
p = @(x) polyval(c,x);
hold on, fplot(p,[-1 1])
%% 
% Surely there must be functions that are more intuitively representative of 
% those points!
% Piecewise Polynomials
% Example 5.1.3

n = 18;
t = linspace(-1,1,n+1)';  
y = t.^2 + t + 0.05*sin(20*t);
clf, plot(t,y,'.')
%% 
% By default |interp1| gives us an interpolating function that is linear between 
% each pair of consecutive nodes.

x = linspace(-1,1,400)';
hold on, plot(x,interp1(t,y,x))
title('Piecewise linear interpolant')    % ignore this line
%% 
% We may instead request a smoother interpolant that is piecewise cubic.

cla
plot(t,y,'.')
plot(x,interp1(t,y,x,'spline'))
title('Piecewise cubic interpolant')   % ignore this line
% Conditioning of Interpolation
% *Example 5.1.4*

figure
n = 18;
t = linspace(-1,1,n+1)';
y = [zeros(9,1);1;zeros(n-9,1)];  % 10th cardinal function
plot(t,y,'.','MarkerSize',15)
title('Data for 10th cardinal function')
%%
hold on
x = linspace(-1,1,400)';
plot(x,interp1(t,y,x,'spline'))
title('Piecewise cubic cardinal function')   % ignore this line
xlabel('x'), ylabel('p(x)')   % ignore this line
hold off
%% 
% The piecewise cubic cardinal function is nowhere greater than one in absolute 
% value. This happens to be true for all the cardinal functions, ensuring a good 
% condition number for the interpolation. But the story for global polynomials 
% is very different.

clf, plot(t,y,'.')    % ignore this line
c = polyfit(t,y,n);
hold on, 
plot(x,polyval(c,x))
title('Polynomial cardinal function')     % ignore this line
xlabel('x'), ylabel('p(x)')   % ignore this line
hold off
%% 
% From the figure we can see that the condition number for polynomial interpolation 
% on these nodes is at least 500. 
% Piecewise Linear Interpolation
% Hat Functions
% *Example 5.2.1*
% 
% Let's define a set of 6 nodes (i.e., $n=5$ in our formulas). 

t = [0 0.075 0.25 0.55 0.7 1]';
%% 
% We plot the hat functions $H_0,\ldots,H_5$. 

for k = 0:5
    fplot(@(x) hatfun(x,t,k),[0 1])
    hold on
end
xlabel('x'), ylabel('H_k(x)'), title('Hat functions')
hold off
% Cardinality Conditions
% *Example 5.2.2*
% 
% We generate a piecewise linear interpolant of $f(x)=e^{\sin 7x}$. 

f = @(x) exp(sin(7*x));
fplot(f,[0 1])
xlabel('x'), ylabel('f(x)')    % ignore this line
title('Function')    % ignore this line
%% 
% First we sample the function to create the data.

t = [0 0.075 0.25 0.55 0.7 1];    % nodes
y = f(t);
hold on, plot(t,y,'.')
title('Function and data')    % ignore this line
%% 
% Now we create a callable function that will evaluate the piecewise linear 
% interpolant at any $x$.

p = plinterp(t,y);
fplot(p,[0 1])
title('Function and PL interpolant')    % ignore this line
%% 
% *Example 5.2.3*
% 
% We measure the convergence rate for piecewise linear interpolation of $e^{\sin 
% 7x}$.

f = @(x) exp(sin(7*x));
x = linspace(0,1,10001)';  % sample the difference at many points
n_ = 2.^(3:10)';
err_ = 0*n_;

for i = 1:length(n_)
    n = n_(i);
    t = linspace(0,1,n+1)';   % interpolation nodes
    p = plinterp(t,f(t));
    err = max(abs( f(x) - p(x) ));
    err_(i) = err;
end
%% 
% Since we expect convergence that is $O(h^2)=O(n^{-2})$, we use a log-log graph 
% of error and expect a straight line of slope $-2$.

figure
loglog( n_, err_, '.-' )
hold on 
loglog( n_, 0.1*(n_/n_(1)).^(-2), '--' )
xlabel('n'), ylabel('||f-p||_\infty')     % ignore this line
title('Convergence of PL interpolation')     % ignore this line
legend('error','2nd order')     % ignore this line
axis tight    % ignore this line
hold off
% Cubic Splines
% *Example 5.3.1*
% 
% For illustration, here is a spline interpolant using just a few nodes. 

f = @(x) exp(sin(7*x));         
fplot(f,[0,1])
t = [0, 0.075, 0.25, 0.55, 0.7, 1]';  % nodes
y = f(t);                            % values at nodes
hold on, 
plot(t,y,'.')
S = spinterp(t,y);
fplot(S,[0,1])
xlabel('x'), ylabel('f(x)')   % ignore this line
title('Function and cubic spline interpolant')   % ignore this line
hold off
%% 
% Now we look at the convergence rate as the number of nodes increases.

x = linspace(0,1,10001)';  % sample the difference at many points
n_ = 2.^(3:0.5:8)';
err_ = 0*n_;

for i = 1:length(n_)
    n = n_(i);
    t = linspace(0,1,n+1)';   % interpolation nodes
    S = spinterp(t,f(t));
    err = norm( f(x) - S(x), inf );
    err_(i) = err;
end
%% 
% Since we expect convergence that is $O(h^4)=O(n^{-4})$, we use a log-log graph 
% of error and expect a straight line of slope $-4$.

clf, loglog( n_, err_, '.-' )
hold on, 
loglog( n_, (n_/n_(1)).^(-4), '--' )
xlabel('n'), ylabel('||f-S||_\infty'), axis tight   % ignore this line
title('Convergence of spline interpolation')   % ignore this line
legend('error','4th order')   % ignore this line
hold off
% Finite Differences
% *Example 5.4.2*
% 
% We try to estimate the derivative of $\cos(x^2)$ at $x=0.5$ using five nodes. 

t = [ 0.35 0.5 0.57 0.6 0.75 ]';   % nodes
f = @(x) cos(x.^2);
dfdx = @(x) -2*x.*sin(x.^2);
exact_value = dfdx(0.5)
%% 
% We have to shift the nodes so that the point of estimation for the derivative 
% is at $x=0$. 

w = fdweights(t-0.5,1);
fd_value = w'*f(t)
%% 
% We can reproduce the weights in the finite difference tables by using equally 
% spaced nodes with $h=1$. For example, here are two one-sided formulas.

format rat
fdweights(0:2,1)
fdweights(-3:0,1)
% Convergence of finite differences
% *Example 5.5.2*
% 
% Let's observe the convergence of the forward-difference formula applied to 
% the function $\sin(e^{x+1})$ at $x=0$. 

f = @(x) sin( exp(x+1) );
FD1 = [ (f(0.1)-f(0))   /0.1
        (f(0.05)-f(0))  /0.05
        (f(0.025)-f(0)) /0.025 ]
%% 
% It's not clear that the sequence is converging. As predicted, however, the 
% errors are cut approximately by a factor of 2 when $h$ is divided by 2.

exact_value = cos(exp(1))*exp(1);
err = exact_value - FD1
%% 
% Asymptotically as $h\to0$, the error is proportional to $h$.
% Higher-order accuracy
% *Example 5.5.4*

f = @(x) sin( exp(x+1) );
exact_value = cos(exp(1))*exp(1);
%% 
% We'll run both formulas in parallel for a sequence of $h$ values.

h = 4.^(-1:-1:-8)';
FD1 = 0*h;  FD2 = 0*h;
for k = 1:length(h)
    FD1(k) = (f(h(k)) - f(0)) / h(k);
    FD2(k) = (f(h(k)) - f(-h(k))) / (2*h(k));
end
%% 
% In each case $h$ is decreased by a factor of 4, so that the error is reduced 
% by a factor of 4 in the first-order method and 16 in the second-order method.

error_FD1 = exact_value-FD1;  
error_FD2 = exact_value-FD2;
table(h,error_FD1,error_FD2)
%% 
% A graphical comparison can be clearer. On a log-log scale, the error should 
% (roughly) be a straight line whose slope is the order of accuracy. However, 
% it's conventional in convergence plots, to show $h$ _decreasing_ from left to 
% right, which negates the slopes.

loglog(h,[abs(error_FD1),abs(error_FD2)],'.-')
hold on, 
loglog(h,[h,h.^2],'--')      % perfect 1st and 2nd order 
set(gca,'xdir','reverse')
xlabel('h'), ylabel('error in f''(0)')     % ignore this line
legend('error FD1','error FD2','order 1','order 2','location','southwest') % ignore this line
title('Convergence of finite differences')   % ignore this line
hold off
% Finite Difference Conditioning
% *Example 5.5.5*
% 
% Let $f(x)=e^{-1.3x}$. We apply finite difference formulas of first, second, 
% and fourth order to estimate $f'(0)=-1.3$.

h = 10.^(-1:-1:-12)';
f = @(x) exp(-1.3*x);
for j = 1:length(h) 
    nodes = h(j)*(-2:2)';
    vals = f(nodes)/h(j);
    fd1(j) = [   0    0 -1    1     0] * vals;
    fd2(j) = [   0 -1/2  0  1/2     0] * vals;
    fd4(j) = [1/12 -2/3  0  2/3 -1/12] * vals;
end
loglog(h,abs(fd1+1.3),'.-'), hold on
loglog(h,abs(fd2+1.3),'.-')
loglog(h,abs(fd4+1.3),'.-')
loglog(h,0.1*eps./h,'--')
set(gca,'xdir','reverse')
legend('FD1','FD2','FD4','O(h^{-1})','location','southeast')  % ignore this line
xlabel('h'), ylabel('total error')     % ignore this line
title('FD error with roundoff'), axis tight    % ignore this line
hold off
%% 
% Again the graph is made so that $h$ decreases from left to right. The errors 
% are dominated at first by truncation error, which decreases most rapidly for 
% the 4th order formula. However, increasing roundoff error eventually equals 
% and then dominates the truncation error as $h$ continues to decrease. As the 
% order of accuracy increases, the crossover point moves to the left (greater 
% efficiency) and down (greater accuracy). 
% Numerical Integration
% *Example 5.6.1*
% 
% The antiderivative of $e^x$ is, of course, itself. That makes evaluation of 
% $\int_0^1 e^x\,dx$ by the Fundamental Theorem trivial.

format long, I = exp(1)-1
%% 
% MATLAB has a built-in |integral| that estimates the value numerically without 
% finding the antiderivative first. As you can see here, it's often just as accurate. 

integral(@(x) exp(x),0,1)
%% 
% The numerical approach is far more robust. For example, $e^{\sin x}$ has no 
% useful antiderivative. But numerically it's no more difficult. 

integral(@(x) exp(sin(x)),0,1)
%% 
% When you look at the graphs of these functions, what's remarkable is that 
% one of these areas is the most basic calculus while the other is almost impenetrable 
% analytically. From a numerical standpoint, they are practically the same problem.

x = linspace(0,1,201)';
subplot(2,1,1), fill([x;1;0],[exp(x);0;0],[1,0.9,0.9])
title('exp(x)')  % ignore this line
ylabel('f(x)')    % ignore this line
subplot(2,1,2), fill([x;1;0],[exp(sin(x));0;0],[1,0.9,0.9])
title('exp(sin(x))')  % ignore this line
xlabel('x'), ylabel('f(x)')    % ignore this line
% Trapezoid Rule
% *Example 5.6.2*
% 
% We approximate the integral of the function $f(x)=e^{\sin 7x}$ over the interval 
% $[0,2]$. 

f = @(x) exp(sin(7*x));
a = 0;  b = 2;
%% 
% In lieu of the exact value, we will use the built-in |integral| function to 
% find an accurate result.

I = integral(f,a,b,'abstol',1e-14,'reltol',1e-14);
fprintf('Integral = %.15f\n',I)
%% 
% Here is the error at $n=40$. 

T = trapezoid(f,a,b,40);
err = I - T
%% 
% In order to check the order of accuracy, we double $n$ a few times and observe 
% how the error decreases.

n = 40*2.^(0:5)';
err = zeros(size(n));
for k = 1:length(n)
    T = trapezoid(f,a,b,n(k));
    err(k) = I - T;
end

table(n,err)
%% 
% Each doubling of $n$ cuts the error by a factor of about 4, which is consistent 
% with second-order convergence. Another check: the slope on a log-log graph should 
% be $-2$.

figure
loglog(n,abs(err),'.-')
hold on, loglog(n,3e-3*(n/n(1)).^(-2),'--')
xlabel('n'), ylabel('error'), axis tight        % ignore this line
title('Convergence of trapezoidal quadrature')  % ignore this line
legend('observed','2nd order')               % ignore this line
hold off
% Extrapolation
% *Example 5.6.3*
% 
% We estimate $\displaystyle\int_0^2 x^2 e^{-2x}\, dx$ using extrapolation. 

f = @(x) x.^2.*exp(-2*x);
a = 0;  b = 2;  format short e
I = integral(f,a,b,'abstol',1e-14,'reltol',1e-14);
%% 
% We start with the trapezoid formula on $n=N$ nodes.

N = 20;       % the coarsest formula
n = N;  h = (b-a)/n;
t = h*(0:n)';   y = f(t);
%% 
% We can now apply weights to get the estimate $T_f(N)$. 

T = h*( sum(y(2:N)) + y(1)/2 + y(n+1)/2 );
err_2nd = I - T
%% 
% Now we double to $n=2N$, but we only need to evaluate $f$ at every other interior 
% node.

n = 2*n;  h = h/2;  t = h*(0:n)';
T(2) = T(1)/2 + h*sum( f(t(2:2:n)) );
err_2nd = I - T
%% 
% As expected for a second-order estimate, the error went down by a factor of 
% about 4. We can repeat the same code to double $n$ again.

n = 2*n;  h = h/2;  t = h*(0:n)';
T(3) = T(2)/2 + h*sum( f(t(2:2:n)) );
err_2nd = I - T
%% 
% Let us now do the first level of extrapolation to get results from Simpson's 
% formula. We combine the elements |T(i)| and |T(i+1)| the same way for $i=1$ 
% and $i=2$.

S = (4*T(2:3) - T(1:2)) / 3;
err_4th = I - S
%% 
% With the two Simpson values $S_f(N)$ and $S_f(2N)$ in hand, we can do one 
% more level of extrapolation to get a 6th-order accurate result.

R = (16*S(2) - S(1)) / 15;
err_6th = I - R
%% 
% If we consider the computational time to be dominated by evaluations of $f$, 
% then we have obtained a result with twice as many accurate digits as the best 
% trapezoid result, at virtually no extra cost.
% Adaptive
% *Example 5.7.1*
% 
% This function gets increasingly oscillatory near the right endpoint.

f = @(x) (x+1).^2.*cos((2*x+1)./(x-4.3));
fplot(f,[0 4],2000)
xlabel('x'), ylabel('f(x)')    % ignore this line
title('A wiggly integrand')
set(gca(),'XScale','linear','YScale','linear')
axis([0,4,-25,25])
%% 
% Accordingly, the trapezoid rule is more accurate on the left half of the interval 
% than on the right half.

n_ = 50*2.^(0:3)';
T_ = [];
for i = 1:length(n_)
    n = n_(i);
    T_(i,1) = trapezoid(f,0,2,n); 
    T_(i,2) = trapezoid(f,2,4,n);
end

left_val = integral(f,0,2,'abstol',1e-14,'reltol',1e-14);
right_val = integral(f,2,4,'abstol',1e-14,'reltol',1e-14);
table(n_,T_(:,1)-left_val,T_(:,2)-right_val,... 
    'variablenames',{'n','left_error','right_error'})
%% 
% Both the picture and the numbers suggest that more nodes should be used on 
% the right half of the interval than on the left half.
%% 
% *Example 5.7.2*

clear
f = @(x) (x+1).^2.*cos((2*x+1)./(x-4.3)); % changed to myfun
I = integral(f,0,4,'abstol',1e-14,'reltol',1e-14)  % 'exact' value
%% 
% We perform the integration and show the nodes selected underneath the curve. 

[Q,t] = intadapt(f,0,4,.001);
fplot(f,[0 4],2000,'k'),
hold on
stem(t,f(t),'.-')
num_nodes = length(t)
title('Adaptive node selection')    % ignore this line
xlabel('x'), ylabel('f(x)')    % ignore this line
hold off
%% 
% The error turns out to be a bit more than we requested. It's only an estimate, 
% not a guarantee.

err = I - Q
%% 
% Let's see how the number of integrand evaluations and the error vary with 
% the requested tolerance.

tol_ = 10.^(-4:-1:-14)';
err_ = 0*tol_;
num_ = 0*tol_;

for i = 1:length(tol_)
    [Q,t] = intadapt(f,0,4,tol_(i));
    err_(i) = I - Q;
    num_(i) = length(t);
end

table(tol_,err_,num_,'variablenames',{'tol','error','f_evals'})
%% 
% As you can see, even though the errors are not less than the estimates, the 
% two columns decrease in tandem. If we consider now the convergence not in $h$ 
% (which is poorly defined) but in the number of nodes actually chosen, we come 
% close to the fourth order accuracy of the underlying Simpson scheme.

clf,  loglog(num_,abs(err_),'.-')
hold on
loglog(num_,0.01*(num_/num_(1)).^(-4),'--')
title('Convergence of adaptive quadrature')   % ignore this line
xlabel('number of nodes'),  ylabel('error')   % ignore this line
legend('error','O(n^{-4})')   % ignore this line
%xlim(num_([1 end]))   % ignore this line
hold off
%% 
% Below are the functions used in this chapter.

function H = hatfun(x,t,k)
% HATFUN   Hat function/piecewise linear basis function.
% Input: 
%   x      evaluation points (vector)
%   t      interpolation nodes (vector, length n+1)
%   k      node index (integer, in 0,...,n)
% Output:
%   H      values of the kth hat function

n = length(t)-1;
k = k+1;  % adjust for starting with index=1

% Fictitious nodes to deal with first, last funcs.
t = [ 2*t(1)-t(2); t(:); 2*t(n+1)-t(n) ];
k = k+1;  % adjust index for the fictitious first node

H1 = (x-t(k-1))/(t(k)-t(k-1));   % upward slope
H2 = (t(k+1)-x)/(t(k+1)-t(k));   % downward slope

H = min(H1,H2);
H = max(0,H);
end

function p = plinterp(t,y)
% PLINTERP   Piecewise linear interpolation.
% Input:
%   t     interpolation nodes (vector, length n+1)
%   y     interpolation values (vector, length n+1)
% Output:
%   p     piecewise linear interpolant (function)

n = length(t)-1;
p = @evaluate;

    % This function evaluates p when called.
    function f = evaluate(x)
        f = 0;
        for k = 0:n
            f = f + y(k+1)*hatfun(x,t,k);
        end
    end

end

function S = spinterp(t,y)
% SPINTERP   Cubic not-a-knot spline interpolation.
% Input:
%   t     interpolation nodes (vector, length n+1)
%   y     interpolation values (vector, length n+1)
% Output:
%   S     not-a-knot cubic spline (function)

t = t(:);  y = y(:);  % ensure column vectors
n = length(t)-1;
h = diff(t);          % differences of all adjacent pairs

% Preliminary definitions.
Z = zeros(n);
I = eye(n);  E = I(1:n-1,:);
J = I - diag(ones(n-1,1),1);
H = diag(h);

% Left endpoint interpolation:
AL = [ I, Z, Z, Z ];
vL = y(1:n);

% Right endpoint interpolation:
AR = [ I, H, H^2, H^3 ];
vR = y(2:n+1);

% Continuity of first derivative:
A1 = E*[ Z, J, 2*H, 3*H^2 ];
v1 = zeros(n-1,1);

% Continuity of second derivative:
A2 = E*[ Z, Z, J, 3*H ];
v2 = zeros(n-1,1);

% Not-a-knot conditions:
nakL = [ zeros(1,3*n), [1,-1, zeros(1,n-2)] ];
nakR = [ zeros(1,3*n), [zeros(1,n-2), 1,-1] ];

% Assemble and solve the full system.
A = [ AL; AR; A1; A2; nakL; nakR ];
v = [ vL; vR; v1; v2; 0 ;0 ];
z = A\v;

% Break the coefficients into separate vectors.
rows = 1:n;
a = z(rows);
b = z(n+rows);  c = z(2*n+rows);  d = z(3*n+rows);
S = @evaulate;

    % This function evaluates the spline when called with a value for x.
    function f = evaulate(x)
        f = zeros(size(x));
        for k = 1:n       % iterate over the pieces
            % Evalaute this piece's cubic at the points inside it.
            index = (x>=t(k)) & (x<=t(k+1));   
            f(index) = polyval( [d(k),c(k),b(k),a(k)], x(index)-t(k) );
        end
    end

end

function w = fdweights(t,m)
%FDWEIGHTS   Fornberg's algorithm for finite difference weights.
% Input:
%   t    nodes (vector, length r+1)
%   m    order of derivative sought at x=0 (integer scalar)
% Output:
%   w    weights for the approximation to the jth derivative (vector)

% This is a compact implementation, not an efficient one. 

r = length(t)-1;
w = zeros(size(t));
for k = 0:r
  w(k+1) = weight(t,m,r,k);
end
end

function c = weight(t,m,r,k)
% Implement a recursion for the weights.
% Input:
%   t   nodes (vector)
%   m   order of derivative sought 
%   r   number of nodes to use from t (<= length(t))
%   k   index of node whose weight is found 
% Output:
%   c   finite difference weight 

if (m<0) || (m>r)        % undefined coeffs must be zero
  c = 0;    
elseif (m==0) && (r==0)  % base case of one-point interpolation
  c = 1;   
else                     % generic recursion 
  if k<r
    c = (t(r+1)*weight(t,m,r-1,k) - ...
        m*weight(t,m-1,r-1,k))/(t(r+1)-t(k+1));
  else
    beta = prod(t(r)-t(1:r-1)) / prod(t(r+1)-t(1:r));
    c = beta*(m*weight(t,m-1,r-1,r-1) - t(r)*weight(t,m,r-1,r-1));
  end
end
end

function [T,t,y] = trapezoid(f,a,b,n)
%TRAPEZOID   Trapezoid formula for numerical integration.
% Input:
%   f     integrand (function)
%   a,b   interval of integration (scalars)
%   n     number of interval divisions
% Output:
%   T     approximation to the integral of f over (a,b)
%   t     vector of nodes used
%   y     vector of function values at nodes

h = (b-a)/n;
t = a + h*(0:n)';
y = f(t);
T = h * ( sum(y(2:n)) + 0.5*(y(1) + y(n+1)) ); 
end

function y=myfun(x)
y=(x+1).^2.*cos((2*x+1)./(x-4.3));
end

function [Q,t] = intadapt(f,a,b,tol)
%INTADAPT   Adaptive integration with error estimation.
% Input:
%   f     integrand (function)
%   a,b   interval of integration (scalars)
%   tol   acceptable error
% Output:
%   Q     approximation to integral(f,a,b)
%   t     vector of nodes used

m = (b+a)/2;
[Q,t] = do_integral(a,f(a),b,f(b),m,f(m),tol);
end

% Use error estimation and recursive bisection. 
function [Q,t] = do_integral(a,fa,b,fb,m,fm,tol)

% These are the two new nodes and their f-values.
xl = (a+m)/2;  fl = myfun(xl);
xr = (m+b)/2;  fr =myfun(xr);
t = [a;xl;m;xr;b];              % all 5 nodes at this level

% Compute the trapezoid values iteratively. 
h = (b-a);
T(1) = h*(fa+fb)/2;
T(2) = T(1)/2 + (h/2)*fm;
T(3) = T(2)/2 + (h/4)*(fl+fr);

S = (4*T(2:3)-T(1:2)) / 3;      % Simpson values
E = (S(2)-S(1)) / 15;           % error estimate
        
if abs(E) < tol*(1+abs(S(2)))   % acceptable error?
    Q = S(2);                   % yes--done
else
    % Error is too large--bisect and recurse. 
    [QL,tL] = do_integral(a,fa,m,fm,xl,fl,tol);
    [QR,tR] = do_integral(m,fm,b,fb,xr,fr,tol);
    Q = QL + QR;
    t = [tL;tR(2:end)];         % merge the nodes w/o duplicate
end        
end