clear
clc
format long
disp ('Método de Simpson 1/3')
syms x
f = input('función: ');
b = input('límite superior: ');
a = input('límite inferior: ');
n = input('número de divisiones (par): ');
h = (b-a)/n;
%integracion normal
i = int(f);
x = b;
r = eval(i);
x = a;
r = r - eval(i);
%metodo de Simpson 1/3
for j=0:n;
if(j == 0)
x = a;
s = eval(f);
elseif(j == n)
x = b;
s = s+eval(f);
elseif(rem(j,2) ==0)
x = x+h;
s = s+2*eval(f);
else
x = x+h;
s = s+4*eval(f);
end
end
s = h*s/3;
fprintf('resultado por simpson 1/3: %f\n',s)
fprintf('resultado por integración normal: %f\n',r)