function [I]=simpson(f,a,b,m)

 

% integral funkcije f v mejah od a do b z m podintervalov

%function [I]=simpson(f,a,b,m)

h=(b-a)/m;

 

x=linspace(a,b,m+1);

y=eval(vectorize(f));

 

I=h/3*(sum(4*y(2:2:m))+sum(2*y(3:2:m))+y(1)+y(m+1));

 

return;