Composite Simpson's Rule

function [outI, outX] = composite_two_simpson (func, a, b, h)
## because of how matlab indexes vectors (1 to "k") instead of (0 to "our n")
## have to change the way we count the starting and ending points of intervals
## h represents distance between mesh points
## can input (b-a)/n in for h if you know how many subintervals you want

if h < (b-a)/2 &&  h > 0

## number of intervals labeled k and initalize other sums
  k = 1+(b-a)/h;
  start_pts = 0;
  end_pts = 0;

## construct outX vector which contains all mesh points
  outX = a:h:b;

## 4f(x(2*i-1)) from 1 to n/2
## instead go from 1 to (k/2)-1
  for i = 1:floor(k/2-1)
    start_pts = start_pts + func(outX(2*i));
    i = i + 1;

  endfor

## 2f(x(2*i)) from 1 to (n/2)-1
## instead go from 2 to (k/2)

  for j = 2:(floor(k/2))
    end_pts = end_pts + func(outX(2*j-1));
    j = j + 1;

  endfor
  outI = h/3*(func(a) + 4*start_pts + 2*end_pts + func(b));

## elseif for simple simpson rule with 2 subintervals
## m is the midpoint
## outI is the approximation of the integral
## outX is the vector containing unique mesh points
elseif h == (b-a)/2
  m = (a+b)/2;
  outI = h/3*(func(a) + 4*func((a+b)/2) + func(b));
  outX = [a, m, b];

else
  print("distance h is too large or a nonpositive distance")

endif

endfunction