Assignment 10 : Romberg Integration

      PROGRAM ROMINT
      DIMENSION R(20,20)
      OPEN(1,'IA10.DAT')
      OPEN(2,'OA10.DAT')
      READ(1,*)A,B,N
      CALL ROMBERG(A,B,N,R)
      END PROGRAM

      SUBROUTINE ROMBERG(A,B,N,R)
      F(X)=SIN(X)
      REAL::R(20,20)
      WRITE(2,*)'ROMBERG INTEGRATION'
      H=B-A
      R(1,1)=(H/2.0)*(F(A)+F(B))
      DO I=2,N
      S=0
      DO K=1,2**(I-2)
      S=S+F(A+((K-0.5)*H))
      END DO
      R(I,1)=0.5*(R(I-1,1)+(H*S))
      H=H/2.0
      END DO
      DO J=2,N
      DO I=J,N
      R(I,J)=R(I,J-1)+((R(I,J-1)-R(I-1,J-1))/(4**(J-1))-1)
      END DO
      END DO
      DO I=1,N
      WRITE(2,20)(R(I,J),J=1,I)
  20  FORMAT(6(2X,F12.6))
      END DO
      END SUBROUTINE