%Numerical Integration methods
%Integrate 1/x^3 from x=a=1 to x=b=3
a=1;
b=3;
n=12;
h=(b-a)/n;
%calculate x and y=f(x) values
xx=a+h*((1:n+1)-1);
yy=1./xx.^3;

%calculate exact solution
exact=(-1/2*(1/b^2-1/a^2));


%method 1 Trapezoidal rule
sum1=h/2*(yy(1)+yy(n+1));
for i=2:n
   sum1=sum1+2*h/2*yy(i);
end


%method 2 simpson's 1/3 rule
sum2=h/3*(yy(1)+yy(n+1));
for i=2:2:n
   sum2=sum2+4*h/3*yy(i);
end
for i=3:2:n-1
   sum2=sum2+2*h/3*yy(i);
end


%method 3 simpson's 3/8 rule
sum3=(3/8*h)*(yy(1)+yy(n+1));
for i=2:3:n-1
    sum3=sum3+3*3/8*h*(yy(i)+yy(i+1));
end
for i=4:3:n-2
    sum3=sum3+2*3/8*h*yy(i);
end

disp(' *Exact soln* *Trapezoidal rule* *Simpson 1/3 rule* *Simpson 3/8 rule*')
[exact sum1 sum2 sum3 ]