function midpt0

%  Integrate  f(x)  over  [a, b]  using composite midpoint rule
%  Compares computed and exact value
%  f(x) given at end of file

a=0; b=1;
exact=(exp(3)-1)/3;

% plot rectangles used for midpoint rule
n=5;
figure(1)
clf
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
% get(gcf)
set(gcf,'Position', [6 1062 627 283])
hold on
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
for in=1:n
    plot([xd(in) xd(in+1)], [f(xd(in)+0.5*h) f(xd(in)+0.5*h)],'--r','LineWidth',1.5)
    plot([xd(in) xd(in)], [0 f(xd(in)+0.5*h)],'--r','LineWidth',1.5)
end
plot([b b], [0 f(xd(n)+0.5*h)],'--r','LineWidth',1.5)
nx=100;
x=linspace(a,b,nx);
for ix=1:nx
    y(ix)=f(x(ix));
end
plot(x,y,'LineWidth',1.5)
box on
xlabel('x-axis')
ylabel('y-axis')
set(gca,'FontSize',16,'FontWeight','bold')

pause

%  calculate integral using composite midpoint rule
nk=6;
interv=[10 20 40 80 160 320];
fprintf('\n Subinervals  Composite Midpoint      Error \n')
for k=1:nk
    n=interv(k);
    xd=linspace(a,b,n+1);
    h=xd(2)-xd(1);
    % calculate I_M
    sum=f(xd(1)+0.5*h);
    for j=2:n
        sum=sum+f(xd(j)+0.5*h);
    end
    I_M(k)=h*sum;
    err(k)=abs(exact-I_M(k));
    fprintf('   n =  %i     I_M = %10.8f     E_M = %8.1e \n',interv(k),I_M(k),err(k));
    pause
end
fprintf('\n')

% plot error
figure(2)
clf
% get(gcf)
set(gcf,'Position', [7 708 625 280])
loglog(interv,err,'--or','MarkerSize',9,'LineWidth',1.5)
hold on
grid on
yticks([1e-5 1e-4 1e-3 1e-2 1e-1])
xlabel('Number of Subintervals (n)')
ylabel('Error')
set(gca,'FontSize',16,'FontWeight','bold')


function y=f(x)
y=exp(3*x);