title: MatLab algorithm for composite Simpson's rule
body: <p>I have tried, just for the fun of it, to write a MatLab-code for the composite Simpson's rule.  As far as I can see, the code is correct, but my answers are not as accurate as I would like.  If I try my code on the function f = cos(x) + e^(x^2), with a = 0, b = 1 and n = 7, my answer is roughly 1,9, when it should be 2,3.  If I use the algorithm available at Wikipedia, I get a very close approximation with n = 7, so my code is obviously not good enough.  If someone can see any mistakes in my code, I would really appreciate it!</p>
  
  <pre><code>function x = compsimp(a,b,n,f)
  % The function implements the composite Simpson's rule
  
  h = (b-a)/n;
  x = zeros(1,n+1);
  x(1) = a;
  x(n+1) = b;
  p = 0;
  q = 0;
  
  % Define the x-vector
  for i = 2:n
      x(i) = a + (i-1)*h;
  end
  
  % Define the terms to be multiplied by 4
  for i = 2:((n+1)/2)
      p = p + (f(x(2*i -2)));
  end
  
  % Define the terms to be multiplied by 2
  for i = 2:((n-1)/2)
      q = q + (f(x(2*i -1)));
  end
  
  % Calculate final output
  x = (h/3)*(f(a) + 2*q + 4*p + f(b));
  </code></pre>