METODO DE SIMPSON 1/3

clc
disp('METODO DE SIMPSON 1/3')
f=inline(input('ingrese la funcion a integrar:','s'));
a=input('ingrese intervalo inferior:');
b=input('ingrese intervalo superior:');
h=input('digite el valor del paso:');
c=0;suma=0;sumav=0;
s=a;
while (s<=b)
    c=c+1;
     s=s+h;
end   
x=a:h:b;
   if(rem(c,2)==0)
     disp('ERROR')
     disp('SE REQUIERE UN NUMERO PAR DE PAREJAS')
     disp('EJECUTE DE NUEVO Y CAMBIE EL PASO')
     break;
 else
   for(i=1:c)
       fi=feval(f,x(i));
       if(i==1 | i==c)
          suma=suma+fi;
      else
          if(rem(i,2)==0)
           sumav=sumav+(fi*4);
       else
           sumav=sumav+(fi*2);
       end 
   end
   end
 end
   suma=suma+sumav; 
   integral=(h/3)*suma;
   t0=feval(f,x(1));
   t1=feval(f,x(2));
   t2=feval(f,x(3));
   t3=feval(f,x(4));
   t4=feval(f,x(5));
   f4=(t4-(4*t3)+(6*t2)-(4*t1)+t0)/(h^4);
   error=-((b-a)/180)*(h^4)*f4;
   integral=integral+error