[Matlab] Algorytm całkowania metodą trapezów

%TrapMethod('x^2 +4',1,6);
%   91.6667

function [ r ] = TrapMethod(f, a ,b )

f=@(x)eval(f);
n = 1000;
s = 0;
dx = (b - a)/n;
i = 1;
    while i < n
        s = s+f(a + i * dx);
        i = i +1;
    end

    s = (s + (f(a) + f(b))/2) * dx;
    r = s;
    disp(s);
end