% Midpoint Rule 
% Input: f - integrant function, a,b - lower and upper boundaries of the
% integration, eps - tolerance
% Output: Approximated integration and number of subintervals (n)
function[fi n]=mid(f,a,b,eps)
s=1; error=1; integ=f((a+b)/2);
while error>eps
    s=s+1; m=2^s; h=(b-a)/m;
    integ1=0;
    for k=1:m/2;
        w=a+(2*k-1)*h;
        integ1=integ1+2*h*f(w);
    end
    error=abs(integ-integ1);
    integ=integ1;
end
fi=integ;
n=2^s;
________________________________________________________________________
clear;
a=0; b=pi/2; eps=10^(-6);
f=inline('(exp(x)-1)/sin(x)');
format long
[integral interval]=mid(f,a,b,eps)
integral =

   3.001595459916656


interval =

        4096