Rewriting trapezoidal to simpson rule in matlab

Question

I'm trying to write a matlab program to calculate an integral by means of trapezoidal and simpsons rule. The program for trapezoidal is as follows:

    function [int, flag, stats] = trapComp(f, a, b, tol, hMin)
        % Initialise variables
        h = b - a;
        n = 1;
        int = h / 2 * (f(a) + f(b));
        flag = 1;

        if nargout == 3
            stats = struct('totalErEst', [], 'totalNrIntervals', [], 'nodesList', []);
        end

        while h > hMin
            h = h / 2;
            n = 2 * n;
            if h < eps % Check if h is not "zero"
                break;
            end

            % Update the integral with the new nodes
            intNew = int / 2;
            for j = 1 : 2 : n
                intNew = intNew + h * f(a + j * h);
            end

            % Estimate the error
            errorEst = 1 / 3 * (int - intNew);
            int = intNew;
            if nargout == 3 % Update stats
                stats.totalErEst = [stats.totalErEst; abs(errorEst)];
                stats.totalNrIntervals = [stats.totalNrIntervals; n / 2];
            end

            if abs(errorEst) < tol
                flag = 0;
                break
            end
        end
    end

Now simpsons rule I cant really quite get around. I know its very similar but I cant seem to figure it out.

This is my simpson code:

function [int, flag, stats] = simpComp(f, a, b, tol, hMin)
    % Initialise variables
    h = b - a;
    n = 1;
    int = h / 3 * (f(a) + 4 * f((a+b)/2) + f(b));
    flag = 1;

    if nargout == 3
        stats = struct('totalErEst', [], 'totalNrIntervals', [], 'nodesList', []);
    end

    while h > hMin
        h = h / 2;
        n = 2 * n;
        if h < eps % Check if h is not "zero"
            break;
        end

        % Update the integral with the new nodes
        intNew = int / 2;
        for j = 1 : 2 : n
            intNew = intNew + h * f(a + j * h);
        end

        % Estimate the error
        errorEst = 1 / 3 * (int - intNew);
        int = intNew;
        if nargout == 3 % Update stats
            stats.totalErEst = [stats.totalErEst; abs(errorEst)];
            stats.totalNrIntervals = [stats.totalNrIntervals; n / 2];
        end

        if abs(errorEst) < tol
            flag = 0;
            break
        end
    end
end

Using this, however, gives an answer for an integral with a larger error than trapezoidal which i feel it shouldnt.

Any help would be appreciated