function fig6_8

%  Integrate  f(x)  over  [a, b]  using midpoint, trap, and simpson's rule
%
%  Inputs:
%        a = left endpoint of interval
%        b = right endpoint of interval
%        n = number of subintervals
%

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

% composite midpoint
intervm=[10 20 40 80 160 320 640];
nkm=length(intervm);
for k=1:nkm
    n=intervm(k);
    M(k)=midpt(a,b,n);
    errm(k)=abs(exact-M(k));
    %    	fprintf(' n =  %i  I_M = %10.8f  E_M = %8.3e \n',interv(k),M(k),err(k));
end

% composite trap
intervt=[10 20 40 80 160 320 640];
nkt=length(intervt);
for k=1:nkt
    n=intervt(k);
    T(k)=trap(a,b,n);
    errt(k)=abs(exact-T(k));
    %        	fprintf(' n =  %i  I_T = %10.8f  E_T = %8.3e \n',interv(k),T(k),err(k));
end

%  composite simpson
intervs=[10 20 40 80 160 320 640];
nks=length(intervs);
for k=1:nks
    n=intervs(k);
    S(k)=simp(a,b,n);
    errs(k)=abs(exact-S(k));
    %        	fprintf(' n =  %i  I_S = %10.8f  E_S = %8.3e \n',interv(k),S(k),err(k));
end

% plot error using loglog plot
clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.06,'P',0.04)
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)

loglog(intervt,errt,'--o','color', [0 0.5 0],'MarkerSize',8,'LineWidth',1.6)
hold on
loglog(intervm,errm,'--db','MarkerSize',8,'LineWidth',1.6)
loglog(intervs,errs,'--*r','MarkerSize',8,'LineWidth',1.6)
grid on
xlabel('n-axis')
ylabel('Error')
legend({' Trap',' Midpt',' Simp'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gca,'/Users/mark/Desktop/mts.eps')


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

%%% composite simpson: n = number of subintervals
function s=simp(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=f(xd(1));
for j=2:2:n
    ff2=f(xd(j+1));
    s=s+4*f(xd(j))+2*ff2;
end
s=h*(s-ff2)/3;

%%% composite trap: n = number of subintervals
function s=trap(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=0.5*f(xd(1));
for j=2:n
    s=s+f(xd(j));
end
s=h*(s+0.5*f(xd(n+1)));

%%% composite midpoint: n = number of subintervals
function s=midpt(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=0;
for j=1:n
    s=s+f(xd(j)+0.5*h);
end
s=h*s;










