program integrate external F real F real poly (4), integral (5) integer size /4/ real xleft, xright real answer1, answer2 logical error call getpoly (poly, size) call printpoly (poly, size) call getintvl (xleft, xright) call integpoly (poly, size, integral) call printpoly (integral, size+1) print *, ' ' answer1 = F(integral,size+1,xright)-F(integral,size+1,xleft) print *, 'exact integration gives ', answer1 call numinteg (poly, size, F, xleft, xright, answer2, error) if (error) then print *, 'cannot integrate numerically!' else print *, 'numerical integration gives ', answer2 endif stop end subroutine getintvl (xleft, xright) real xleft, xright print *, 'leftmost point of interval of integration?' read *, xleft print *, 'rightmost point of interval of integration?' read *, xright return end