function trap0

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

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

% plot trapezoids used for trapezoidal 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)) f(xd(in+1))],'--r','LineWidth',1.4)
    plot([xd(in) xd(in)], [0 f(xd(in))],'--r','LineWidth',1.4)
end
plot([b b], [0 f(xd(n+1))],'--r','LineWidth',1.4)
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 trapezoidal rule
nk=6;
interv=[10 20 40 80 160 320];
fprintf('\n Subinervals  Composite Trapezoidal      Error \n')
for k=1:nk
    n=interv(k);
    % calculate I_T
    I_T(k)=trap(a,b,n);
    err(k)=abs(exact-I_T(k));
    fprintf('   n =  %i     I_M = %10.8f      E_M = %8.1e \n',interv(k),I_T(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
% include error using composite midpoint
for k=1:nk
    n=interv(k);
    I_M(k)=midpt(a,b,n);
    err_M(k)=abs(exact-I_M(k));
end
loglog(interv,err_M,'--ob','MarkerSize',9,'LineWidth',1.5)
legend({' Trap',' Midpoint'},'Location','NorthEast','FontSize',16,'FontWeight','bold')
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 g=midpt(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
sum=f(xd(1)+0.5*h);
for j=2:n
    sum=sum+f(xd(j)+0.5*h);
end
g=h*sum;

function g=trap(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
sum=0.5*f(xd(1));
for j=2:n
    sum=sum+f(xd(j));
end
g=h*(sum+0.5*f(xd(n+1)));



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