function integrateS

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

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

ic=0;
for k=2:10
    % n =  number of subintervals
    n=2^k;
    if k==10
        n=1000;
    end
    ic=ic+1;

    % midpoint
    evalsM(ic)=n;
    I_M(ic)=midpt(a,b,n);
    errM(ic)=abs(exact-I_M(ic));
    
    % simpson
    evalsS(ic)=n+1;
    I_S(ic)=simp(a,b,n);
    errS(ic)=abs(exact-I_S(ic));
    
    %  2-point gaussian
    evalsG(ic)=2*n;
    I_G(ic)=gauss2(a,b,n);
    errG(ic)=abs(exact-I_G(ic));
    
end

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])
loglog(evalsM,errM,'--ob','MarkerSize',8,'LineWidth',1.3)
hold on
loglog(evalsS,errS,'-*r','MarkerSize',9,'LineWidth',1.3)
loglog(evalsG,errG,'-.+k','MarkerSize',9,'LineWidth',1.3)
grid on
xlabel('Number of f(x) Evaluations')
ylabel('Error')
axis([1 1.01e3 1e-16 1]);
set(gca,'ytick',[1e-16 1e-11 1e-6 1e-1])
set(gca,'YMinorGrid','off')
legend({' Midpoint',' Simpson',' Gauss'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')



%  composite 2-point gaussian
function g=gauss2(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
c1=0.5*h*(1-1/sqrt(3));
c2=0.5*h*(1+1/sqrt(3));
sum= 0;
for j=1:n
    sum=sum+f(xd(j)+c1)+f(xd(j)+c2);
end
g=0.5*h*sum;

%  composite Simpson
function g=simp(a,b,n)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
sum=f(a);
for j=2:2:n
    ff2=f(xd(j+1));
    sum=sum+4*f(xd(j))+2*ff2;
end
g=h*(sum-ff2)/3;

%  composite midpoint
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 y=f(x)
y=exp(3*x);