Simpson's 3/8 Rule

Could someone please take a look at my code. When I execute this code it says: ??? Index exceeds matrix dimensions. Plus it gave me the wrong value. I don't know what went wrong. Thank you in advance.

% Inputs
%       This is the only place in the program where the user makes the changes
%       based on their wishes

% f(x), the function to integrate

f= @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400*x.^5 ;

% a, the lower limit of integration

a= 0 ;

% b, the upper limit of integration

b= 0.8 ;

% n, the maximum number of segments

n=5 ;

%**********************************************************************

% Displays what inputs are being used

disp(sprintf('\n\n********************************Input Data**********************************\n'))
disp(sprintf('     f(x), integrand function'))
disp(sprintf('     a = %g, lower limit of integration ',a))
disp(sprintf('     b = %g, upper limit of integration ',b))
disp(sprintf('     n = %g, number of subdivision',n))
format short g

% Determine the exact solution

exact = quad(f,a,b) ;

% Determine number of different segments to consider

nstep = floor(log2(n)) ;

for i=0:nstep-1

      % Determine number of segments and h
      NN(i+1)=2^(i+1) ;
      h=(b-a)/NN(i+1) ;

      % Use the rule to find the approximation
      integral = f(a) + f(b) ;
      for j=1:2:NN(i+1)-1
          integral = integral + 3*f(a+h*j) ;
      end
      for j=2:2:NN(i+1)-2
          integral = integral + 3*f(a+h*j) ;
      end
      for j=3:2:NN(i+1)-3
          integral = integral + 2*f(a+h*j) ;
      end
      integral = integral * 3 * h/8 ;
      YY(i+1)=integral ;

      % Compute Errors
      Et(i+1)=exact-integral ;
      Etabs(i+1)=abs((integral-exact)/exact) ;
      if(i > 0)
          Ea(i+1)=YY(i+1)-YY(i) ;
          Eaabs(i+1)=abs((YY(i+1)-YY(i))/YY(i)) ;
          SD(i+1)=floor((2-log10(Eaabs(i+1)/0.5))) ;
          if(SD(i+1)<0)
              SD(i+1)=0
          end
      else
          Ea(1)=0 ;
          Eaabs(1)=0 ;
          SD(1)=0 ;
      end

end

% The following builds and displays a table of values

disp(sprintf('\n\n************************Table of Values******************************\n'))

disp('        Approx        True     Relative      Approx   Rel Appr   Sig   ')
disp(' n     Integral       Error   True Error     Error     Error   Digits ')
disp('-----------------------------------------------------------------------')

for i=1:nstep

      if(i > 1)
          if(exact || YY(i) > 0)
              disp(sprintf('%4i  %+1.3e  %+1.3e  %+1.3e  %+1.3e  %+1.3e  %2i',NN(i),YY(i),Et(i),Etabs(i),Ea(i),Eaabs(i),SD(i) ))
          else
              disp(sprintf('%4i  %+1.3e  %+1.3e     n/a      %+1.3e      n/a     n/a',NN(i),YY(i),Etabs(i),Ea(i)))
          end
      else
          disp(sprintf('%4i  %+1.3e  %+1.3e  %+1.3e      n/a         n/a     n/a',NN(i),YY(i),Et(i),Etabs(i)))
      end

end
disp('-----------------------------------------------------------------------')

% This function detects information about your
% screensize and tries to then place/size the graphs accordingly.
scnsize = get(0,'ScreenSize');

% Graph: The following code sets up the graphical depiction of the method.
x(1)=a ;
y(1)=f(a) ;
hold on
for i=1:n
    x(i+1)=a+i*h ;
    y(i+1)=f(x(i+1)) ;
end

for i=1:2:n
    bg = i ;
    ed = i + 2 ;
    coeffs = polyfit(x(bg:ed), y(bg:ed), 2);
    poly_x1 = [x(bg):(x(ed) - x(bg))/1000:x(ed)];
    poly_y1 = coeffs(1)*poly_x1.^2 + coeffs(2)*poly_x1 + coeffs(3);
    poly_x1 = [poly_x1(1) poly_x1 poly_x1(end)];
    poly_y1 = [0 poly_y1 0];
    fill(poly_x1, poly_y1, 'y')
end
xrange=a:(b-a)/1000:b;
plot(xrange,f(xrange),'k','Linewidth',2)
title('Integrand function and Graphical Depiction of Simpson''s 3/8 Rule')

Answer 1

[Oleg Komarov]

The problem is in the last LOOP, you set at the last iteration:

ed = i + 2;

When i = 5, ed = 7 but x has only 6 elements.

EDIT Also note that on FEX just recently somebody submitted: http://www.mathworks.com/matlabcentral/fileexchange/33493-simpsons-13-and-38-rules

[Akmal]

But when I changed it to ed = i + 1, I got this error :

Warning: Polynomial is not unique; degree >= number of data points.

> In polyfit at 71

Warning: Polynomial is not unique; degree >= number of data points.

> In polyfit at 71

Warning: Polynomial is not unique; degree >= number of data points.

> In polyfit at 71

I'm at my wit's end!

[Oleg Komarov]

Try:

for i=2:2:n-1

bg = i-1 ;

ed = i + 1 ;

coeffs = polyfit(x(bg:ed), y(bg:ed), 2);

...

[Akmal]

It works! Thank you so much!