function Imp = midpoint(a,b,M,f)
% USE: midpoint(0,1,10,'cos')
% or
% midpoint(0,1,10,@fun)
H = (b-a)/M;
x = linspace(a+H/2,b-H/2,M);
fmp = feval(f,x);
Imp=H*sum(fmp);