function fig6_9

%  Integrate  f(x)  over  [a, b]  using composite simpson and simpson-romberg

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

%  interv = number of subinternals
interv=[4 8 16 32 64 128];
nk=length(interv);
tol=1e-10;
n=2;
intervals(1)=n;
S(1)=simp(a,b,2);
errS(1)=abs(S(1)-exact)/exact;
err=10*tol;
R(1)=1;  errR(1)=-1;  % this is a fake value
k=1;
while err>tol
    k=k+1;
    n=2*n;
    intervals(k)=n;
    S(k)=simp(a,b,n);
    R(k)=(16*S(k)-S(k-1))/15;
    err=abs(1-R(k)/R(k-1));
    errR(k)=abs(R(k)-exact)/exact;
    errS(k)=abs(S(k)-exact)/exact;
end

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(intervals,errS,'--ob','MarkerSize',8,'LineWidth',1.6)
hold on
loglog(intervals(2:k),errR(2:k),'--*r','MarkerSize',8,'LineWidth',1.6)
grid on
xlabel('n-axis')
ylabel('Error')
legend({' Simp',' Simp-Romb'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
axis([1 2e2 1e-14 1e-1])
set(gca,'ytick',[1e-14 1e-10 1e-6 1e-2])
set(gca,'YMinorGrid','off')
set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gca,'/Users/mark/Desktop/rom2.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;














