% FIGURE 1.6.2
      function YOUT = DSPLN(X,Y,N,YXX1,YXXN,XOUT,NOUT)
%
%  FUNCTION DSPLN FITS AN INTERPOLATORY CUBIC SPLINE THROUGH THE
%    POINTS (X(I),Y(I)), I=1,...,N, WITH SPECIFIED SECOND DERIVATIVES
%    AT THE END POINTS, AND EVALUATES THIS SPLINE AT THE OUTPUT POINTS
%    XOUT(1),...,XOUT(NOUT).
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    X      - A VECTOR OF LENGTH N CONTAINING
%             THE X-COORDINATES OF THE DATA
%             POINTS.
%
%    Y      - A VECTOR OF LENGTH N CONTAINING
%             THE Y-COORDINATES OF THE DATA
%             POINTS.
%
%    N      - THE NUMBER OF DATA POINTS
%             (N >= 3).
%
%    YXX1   - THE SECOND DERIVATIVE OF THE
%             CUBIC SPLINE AT X(1).
%
%    YXXN   - THE SECOND DERIVATIVE OF THE
%             CUBIC SPLINE AT X(N).  (YXX1=0
%             AND YXXN=0 GIVES A NATURAL
%             CUBIC SPLINE)
%
%    XOUT   - A VECTOR OF LENGTH NOUT CONTAINING
%             THE X-COORDINATES AT WHICH THE
%             CUBIC SPLINE IS EVALUATED.  THE
%             ELEMENTS OF XOUT MUST BE IN
%             ASCENDING ORDER.
%
%    YOUT   -                                   A VECTOR OF LENGTH NOUT.
%                                               YOUT(I) CONTAINS THE
%                                               VALUE OF THE SPLINE
%                                               AT XOUT(I).
%
%    NOUT   - THE NUMBER OF OUTPUT POINTS.
%
%-----------------------------------------------------------------------
      SIG(1) = YXX1;
      SIG(N) = YXXN;
%                             SET UP TRIDIAGONAL SYSTEM SATISFIED
%                             BY SECOND DERIVATIVES (SIG(I)=SECOND
%                             DERIVATIVE AT X(I)).
      for I=1:N-2
         HI   = X(I+1)-X(I);
         HIP1 = X(I+2)-X(I+1);
         R(I) = (Y(I+2)-Y(I+1))/HIP1 - (Y(I+1)-Y(I))/HI;
         A(I,1) = HI/6.0;
         A(I,2) = (HI + HIP1)/3.0;
         A(I,3) = HIP1/6.0;
         if (I == 1)
            R(1)   = R(1)   - HI/  6.0*SIG(1);
         end
         if (I == N-2)
            R(N-2) = R(N-2) - HIP1/6.0*SIG(N);
         end
      end
%                             CALL DBAND TO SOLVE TRIDIAGONAL SYSTEM
      NLD = 1;
      NUD = 1;
      SIG(2:N-1) = DBAND(A,N-2,NLD,NUD,R);
%                             CALCULATE COEFFICIENTS OF CUBIC SPLINE
%                             IN EACH SUBINTERVAL
      for I=1:N-1
         HI = X(I+1)-X(I);
         COEFF(1,I) = Y(I);
         COEFF(2,I) = (Y(I+1)-Y(I))/HI - HI/6.0*(2*SIG(I)+SIG(I+1));
         COEFF(3,I) = SIG(I)/2.0;
         COEFF(4,I) = (SIG(I+1)-SIG(I))/(6.0*HI);
      end
      L = 1;
      for I=1:NOUT
%                             FIND FIRST VALUE OF J FOR WHICH X(J+1) IS
%                             GREATER THAN OR EQUAL TO XOUT(I).  SINCE
%                             ELEMENTS OF XOUT ARE IN ASCENDING ORDER,
%                             WE ONLY NEED CHECK THE KNOTS X(L+1)...X(N)
%                             WHICH ARE GREATER THAN OR EQUAL TO
%                             XOUT(I-1).
         for J=L:N-1
            JSAVE = J;
            if (X(J+1) >= XOUT(I))
               break
            end
         end
         L = JSAVE;
%                             EVALUATE CUBIC SPLINE IN INTERVAL
%                             (X(L),X(L+1))
         P = XOUT(I)-X(L);
         YOUT(I) = COEFF(1,L)     + COEFF(2,L)*P       ...
                 + COEFF(3,L)*P*P + COEFF(4,L)*P*P*P;
      end