[Matlab] Algorytm całkowania metodą Simpsona

%>> Tested('x^2 +4',1,6,1)
%ans =   91.6667

function [ o ] = Tested( f, a, b, d)
syms x;
f = @(x)eval(f);
x = [a:d:b];
s = 0;
%
    for i = 1:(length(x) - 1)
       p = ( f(x(i)) + 4 * f((x(i)+ d/2)) + f(x(i)+d) ) * d / 6;
       s = s + p;
    end
o = s;
end