% simpson Composite Simpson's rule
%
% Output: int = approximate value of the integral from a to b of f(x)*dx

n = 5;     % input total number of nodes must be odd number 5 7 9 ....
a = 8;     % input lower limit
b = 30;    % input upper limit
h = (b-a)/(n-1);     %  stepsize
x = a:h:b;           %  divide the interval
fun = 2000*ln(140000./(140000-2100*x))-9.8*x; % input function


f = eval("fun");    %  evaluate integrand

% output
int = (h/3)*( f(1) + 4*sum(f(2:2:n-1)) + 2*sum(f(3:2:n-2)) + f(n) );
%            f(a)       f_even              f_odd            f(b)
write("integral = ",int);