function simp0

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

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

% plot piecewise quadratic approximations
n=2;
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:2:n-1
    % p_2(x)
    ns=10;
    xs=linspace(xd(in),xd(in+2),ns);
    for ix=1:ns
        s(ix)=p2(xs(ix),xd(in),xd(in+1),xd(in+2),f(xd(in)),f(xd(in+1)),f(xd(in+2)));
    end   
    plot(xs,s,'--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 Simpson 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_S(k)=simp(a,b,n);
    err(k)=abs(exact-I_S(k));
    fprintf('   n =  %i     I_M = %10.8f      E_M = %8.1e \n',interv(k),I_S(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 midpoint and trap
for k=1:nk
    n=interv(k);
    I_M(k)=midpt(a,b,n);
    err_M(k)=abs(exact-I_M(k));
    I_T(k)=trap(a,b,n);
    err_T(k)=abs(exact-I_T(k));
end
loglog(interv,err_T,'--ob','MarkerSize',9,'LineWidth',1.5)
loglog(interv,err_M,'--om','MarkerSize',9,'LineWidth',1.5)
legend({' Simpson',' 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 y=p2(x,x1,x2,x3,y1,y2,y3)
l1=(x-x2)*(x-x3)/( (x1-x2)*(x1-x3) );
l2=(x-x1)*(x-x3)/( (x2-x1)*(x2-x3) );
l3=(x-x1)*(x-x2)/( (x3-x1)*(x3-x2) );
y=y1*l1+y2*l2+y3*l3;

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;

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);