function ypp = spline_cubic_set ( n, t, y, ibcbeg, ybcbeg, ibcend, ybcend )

%*****************************************************************************80
%
%% SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline.
%
%  Discussion:
%
%    For data interpolation, the user must call SPLINE_CUBIC_SET to
%    determine the second derivative data, passing in the data to be
%    interpolated, and the desired boundary conditions.
%
%    The data to be interpolated, plus the SPLINE_CUBIC_SET output,
%    defines the spline.  The user may then call SPLINE_CUBIC_VAL to
%    evaluate the spline at any point.
%
%    The cubic spline is a piecewise cubic polynomial.  The intervals
%    are determined by the "knots" or abscissas of the data to be
%    interpolated.  The cubic spline has continous first and second
%    derivatives over the entire interval of interpolation.
%
%    For any point T in the interval T(IVAL), T(IVAL+1), the form of
%    the spline is
%
%      SPL(T) = A(IVAL)
%             + B(IVAL) * ( T - T(IVAL) )
%             + C(IVAL) * ( T - T(IVAL) )^2
%             + D(IVAL) * ( T - T(IVAL) )^3
%
%    If we assume that we know the values Y(*) and YPP(*), which represent
%    the values and second derivatives of the spline at each knot, then
%    the coefficients can be computed as:
%
%      A(IVAL) = Y(IVAL)
%      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
%        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
%      C(IVAL) = YPP(IVAL) / 2
%      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
%
%    Since the first derivative of the spline is
%
%      SPL'(T) =     B(IVAL)
%              + 2 * C(IVAL) * ( T - T(IVAL) )
%              + 3 * D(IVAL) * ( T - T(IVAL) )^2,
%
%    the requirement that the first derivative be continuous at interior
%    knot I results in a total of N-2 equations, of the form:
%
%      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
%      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL)
%
%    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
%
%      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
%      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
%      + YPP(IVAL-1) * H(IVAL-1)
%      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
%      =
%      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
%      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
%
%    or
%
%      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
%      + YPP(IVAL) * H(IVAL)
%      =
%        6 * ( Y(IVAL+1) - Y(IVAL)   ) / H(IVAL)
%      - 6 * ( Y(IVAL)   - Y(IVAL-1) ) / H(IVAL-1)
%
%    Boundary conditions must be applied at the first and last knots.
%    The resulting tridiagonal system can be solved for the YPP values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663.
%
%  Parameters:
%
%    Input, integer N, the number of data points; 
%    N must be at least 2.
%
%    Input, real T(N), the points where data is specified.
%    The values should be distinct, and increasing.
%
%    Input, real Y(N), the data values to be interpolated.
%
%    Input, integer IBCBEG, the left boundary condition flag:
%    0: the spline should be a quadratic over the first interval;
%    1: the first derivative at the left endpoint should be YBCBEG;
%    2: the second derivative at the left endpoint should be YBCBEG;
%    3: Not-a-knot: the third derivative is continuous at T(2).
%
%    Input, real YBCBEG, the left boundary value, if needed.
%
%    Input, integer IBCEND, the right boundary condition flag:
%    0: the spline should be a quadratic over the last interval;
%    1: the first derivative at the right endpoint should be YBCEND;
%    2: the second derivative at the right endpoint should be YBCEND;
%    3: Not-a-knot: the third derivative is continuous at T(N-1).
%
%    Input, real YBCEND, the right boundary value, if needed.
%
%    Output, real YPP(N), the second derivatives of the cubic spline.
%

%
%  Check.
%
  if ( n <= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_CUBIC_SET - Fatal error!\n' );
    fprintf ( 1, '  The number of knots must be at least 2.\n' );
    fprintf ( 1, '  The input value of N = %d\n', n );
    error ( 'SPLINE_CUBIC_SET - Fatal error!' );
  end

  for i = 1 : n - 1
    if ( t(i+1) <= t(i) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPLINE_CUBIC_SET - Fatal error!\n' );
      fprintf ( 1, '  The knots must be strictly increasing, but\n' );
      fprintf ( 1, '  T(%d) = %f\n', i,   t(i)   );
      fprintf ( 1, '  T(%d) = %f\n', i+1, t(i+1) );
      error ( 'SPLINE_CUBIC_SET - Fatal error!' );
    end
  end
%
%  Set space for the sparse matrix A.
%
  a = sparse ( [], [], [], n, n, 3 * n );
  b = zeros ( n, 1 );
%
%  Set the first equation.
%
  if ( ibcbeg == 0 )
    b(1) = 0.0;
    a(1,1) = 1.0;
    a(1,2) = -1.0;
  elseif ( ibcbeg == 1 )
    b(1) = ( y(2) - y(1) ) / ( t(2) - t(1) ) - ybcbeg;
    a(1,1) = ( t(2) - t(1) ) / 3.0;
    a(1,2) = ( t(2) - t(1) ) / 6.0;
  elseif ( ibcbeg == 2 )
    b(1) = ybcbeg;
    a(1,1) = 1.0;
    a(1,2) = 0.0;
  elseif ( ibcbeg == 3 )
    b(1) = 0.0;
    a(1,1) = - ( t(3) - t(2) );
    a(1,2) =   ( t(3)        - t(1) );
    a(1,3) = - (        t(2) - t(1) );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_CUBIC_SET - Fatal error!\n' );
    fprintf ( 1, '  The boundary flag IBCBEG must be 0, 1 or 2.\n' );
    fprintf ( 1, '  The input value is IBCBEG = %d\n', ibcbeg );
    error ( 'SPLINE_CUBIC_SET - Fatal error!' );
  end
%
%  Set the intermediate equations.
%
  for i = 2 : n - 1
    b(i) = ( y(i+1) - y(i) ) / ( t(i+1) - t(i) ) ...
           - ( y(i) - y(i-1) ) / ( t(i) - t(i-1) );
    a(i,i-1) = ( t(i) - t(i-1) ) / 6.0;
    a(i,i) = ( t(i+1) - t(i-1) ) / 3.0;
    a(i,i+1) = ( t(i+1) - t(i) ) / 6.0;
  end
%
%  Set the last equation.
%
  if ( ibcend == 0 )
    b(n) = 0.0;
    a(n,n-1) = -1.0;
    a(n,n) = 1.0;
  elseif ( ibcend == 1 )
    b(n) = ybcend - ( y(n) - y(n-1) ) / ( t(n) - t(n-1) );
    a(n,n-1) = ( t(n) - t(n-1) ) / 6.0;
    a(n,n) = ( t(n) - t(n-1) ) / 3.0;
  elseif ( ibcend == 2 )
    b(n) = ybcend;
    a(n,n-1) = 0.0;
    a(n,n) = 1.0;
  elseif ( ibcend == 3 )
    b(n) = 0.0;
    a(n,n-2) = - ( t(n) - t(n-1) );
    a(n,n-1) =   ( t(n)          - t(n-2) );
    a(n,n) =   - (        t(n-1) - t(n-2) );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_CUBIC_SET - Fatal error!\n' );
    fprintf ( 1, '  The boundary flag IBCEND must be 0, 1 or 2.\n' );
    fprintf ( 1, '  The input value is IBCEND = %d\n', ibcend );
    error ( 'SPLINE_CUBIC_SET - Fatal error!' );
  end
%
%  Special case:
%    N = 2, IBCBEG = IBCEND = 0.
%
  if ( n == 2 && ibcbeg == 0 && ibcend == 0 )

    ypp(1) = 0.0;
    ypp(2) = 0.0;
%
%  Solve the sparse linear system.
%
  else

    ypp = a \ b;

  end

  return
end
