  function PiViaSpline()
% Compute a spline interpolant of sin(x) across [0,2pi].
% Its second derivative will look like -sin(x).
% Its square will look like sin(x)^2 whose integral form 0 to 2pi is about 
% equal to pi.

xVals = linspace(0,2*pi,400)';
for n=[10 20 50 100 500 1000 2000]
    x = linspace(0,2*pi,n)';
    y = sin(x);
    % Compute a spline s(x) interpolant of sin(x) on [0,2pi]
    S = spline(x,y);
    [x,rho] = unmkpp(S);
    % Using the notation xk = x(k), the kth cubic is
    %       rho(k,1)(x-xk)^3 + rho(k,2)(x-xk)^2 + rho(k,3)(x-xk) + rho(k,4) 
    %
    % Its second derivative is
    %       6*rho(k,1)(x-xk) + 2*rho(k,2)
    %
    % Its square is
    %      36*rho(k,1)(x-xk)^2 + 24*rho(k,1)*rho(k,2) + 4*rho(k,2)^2
    %
    % Compute the representation of [s''(x)]^2 and display
    rho2 = [36*rho(:,1).^2   24*rho(:,1).*rho(:,2)  4*rho(:,2).^2];
    ddS = mkpp(x,rho2);
    xVals = linspace(0,2*pi,400)';
    plot(xVals,ppval(ddS,xVals))
    
    % Integrate [s"(x)]^2 from x(1) to x(n)
    % The k-th polynomial is the quadratic
    %       rho2(k,1)*(x-xk)^2 + rho2(k,2)(x-xk) + rho2(k,3)
    % Its integral from x(k) to x(k+1) is
    %       (1/6)rho2(k,1)*delx^3 + (1/2)rho2(k,2)*delx^2 + rho2(k,3)*delx
    % where delx = x(k+1)-x(k)
    delx = x(2:n)'-x(1:n-1)';
    vec_of_integrals = (((rho2(:,1)/6).*delx) + (rho2(:,2)/2).*delx + rho2(:,3)).*delx;
    int_of_ddS = sum(vec_of_integrals);
    title(sprintf('n = %3d  integral = %8.4f',n,int_of_ddS),'fontsize',14)
    shg
    pause
end
    
    