program LabSpline
        integer :: n
        real :: res, arg
        real, dimension(7) :: x, y, b, c, d

        n = 7
        h = 2/25.0
        x = (/-4, -2, 0, 2, 4, 6, 8/) 
        y = (/0, -2, -1, 4, 5, 2, 0/)
           
         call SPLINE(n, x, y, b, c, d)
         
         do i = 1, 151, 1
           arg = -4 + (i-1) * h
           res = SEVAL(n, arg, x, y, b, c, d)
           write(*,*) arg, res
         end do

     end program LabSpline

      SUBROUTINE SPLINE(N,X,Y,B,C,D)

      INTEGER N

      REAL X(N),Y(N),B(N),C(N),D(N)

      INTEGER NM1,IB,I

      REAL T

      NM1=N-1

      IF(N.LT.2) RETURN

      IF(N.LT.3) GO TO 50


      D(1)=X(2)-X(1)

      C(2)=(Y(2)-Y(1))/D(1)

      DO 10 I=2,NM1

      D(I)=X(I+1)-X(I)

      B(I)=2.*(D(I-1)+D(I))

      C(I+1)=(Y(I+1)-Y(I))/D(I)

      C(I)=C(I+1)-C(I)

      10 CONTINUE


      B(1)=-D(1)

      B(N)=-D(N-1)

      C(1)=0.

      C(N)=0.

      IF(N.EQ.3) GO TO 15

      C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))

      C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))

      C(1)=C(1)*D(1)**2/(X(4)-X(1))

      C(N)=-C(N)*D(N-1)**2/(X(N)-X(N-3))


      15 DO 20 I=2,N

      T=D(I-1)/B(I-1)

      B(I)=B(I)-T*D(I-1)

      C(I)=C(I)-T*C(I-1)

      20 CONTINUE

      C(N)=C(N)/B(N)

      DO 30 IB=1,NM1

      I=N-IB

      C(I)=(C(I)-D(I)*C(I+1))/B(I)

      30 CONTINUE

      B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.*C(N))

      DO 40 I=1,NM1

      B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.*C(I))

      D(I)=(C(I+1)-C(I))/D(I)

      C(I)=3.*C(I)

      40 CONTINUE

      C(N)=3.*C(N)

      D(N)=D(N-1)

      RETURN

      50 B(1)=(Y(2)-Y(1))/(X(2)-X(1))

      C(1)=0.

      D(1)=0.

      B(2)=B(1)

      C(2)=0.

      D(2)=0.

      RETURN

      END

      REAL FUNCTION SEVAL(N,U,X,Y,B,C,D)

      INTEGER N

      REAL U,X(N),Y(N),B(N),C(N),D(N)

      INTEGER I,J,K

      REAL DX

      DATA I/1/

      IF(I.GE.N) I=1

      IF(U.LT.X(I)) GO TO 10

      IF(U.LE.X(I+1)) GO TO 30

      10 I=1

      J=N+1

      20 K=(I+J)/2

      IF(U.LT.X(K))J=K

      IF(U.GE.X(K))I=K

      IF(J.GT.I+1)GO TO 20

      30 DX=U-X(I)

      SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))

      RETURN

      END