PROGRAM prog3B
!Initialise the program and give it a name

	IMPLICIT NONE
	!Assume nothing about the variables used within this program
	INTEGER :: n, i
	!Declare integer vars for the number of subintervals and a counting var
	REAL(16) :: h, sum1, sum2, trapz, int, error
	!Delcare high-precision real vars for the step size, sums of terms, trapezium rule result, integration
	!result and error between the two. High-precision is used as the deviation and step size may be very small
	
	PRINT*, 'How many subintervals would you like to calculate the integral over?'
	!Prompt the user to input the number of subintervals
	READ*, n
	!Read the value into variable n
	
	h = 1.0_16/(1.0_16*n)
	!Calculate the step size as interval width / number of subintervals
	sum1 = 3.0_16*h/4.0_16
	!Calculate the (h/2)*(f(x0)+f(xn)) where we're integrating over [x0, xn]
	sum2 = 0.0_16
	!Set the initial value of the sum of intermediate terms to be 0 as would be the case if n = 1
	
	DO i = 1, n-1
	!Open a definite DO loop to sum the intermediate terms running from 1 to n-1 (we don't want x0 and xn)
		sum2 = sum2 + h/(1.0_16+(i*h)**2)
		!Calculate the sum of all previous intermediate terms plus the next
		![sum of previous terms + (h/2)*2*f(xi)]
	END DO
	!End the DO loop
	
	trapz = sum1 + sum2
	!Calculate the result of the trapezium rule with step size h
	int = ATAN(1.0_16)
	!Calculate the value of the integral of 1/(1+x^2) over x = [0,1]
	error = int - trapz
	!Calculate the deviation of the trapezium rule result from the integration result
	
	PRINT '(a38, e11.5, a47)', 'By trapezium rule with step size, h = ', h, &
		', the integral of 1/(1+x^2) between 0 and 1 is:'
	PRINT '(f12.10)', trapz
	!Print the step size and trapezium rule result with limited output precision

	PRINT*, 'The integral of 1/(1+x^2) between 0 and 1 is:'
	PRINT '(f12.10)', int
	!Print out the actual integration result with limited ouput precision

	PRINT*, 'The error between integration of the function and the trapezium rule is:'
	PRINT '(e16.10)', error
	!Print out the error with limited output precision
	
END PROGRAM prog3B
!End the program
