Are you looking for something like this?
a = 1; b = 3;
f = @(x) exp(x);
n = 6;
Ysimp = simp13(f,a,b,n)
Yexact = exp(b)-exp(a)
function [duck] = simp13(func, a, b, n)
%(b-a)/6 * f(x0) + 4f(x1) + f(x2)
if b <= a; error("b cannot be greater than a"); end
% n equals number of times rule is implimented
l = b-a;
c = l/(2*n);
xvals = a:c:b;
duck = 0;
count = 1;
for i = 1:n
duck = duck + (xvals(:,count + 2) - xvals(:,count))/6*(func(xvals(:,count)) + 4*func(xvals(:,count + 1)) + func(xvals(:,count + 2)));
count = count + 2;
disp([count, duck])
end
end
