METODO DE SIMPSON 3/8

clc
disp('METODO DE SIMPSON 3/8')
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-1),3)~=0)
     disp('nnERROR')
     disp('SE REQUIERE UN NUMER0 DE PAREJAS DIVISIBLE POR  3')
     disp('EJECUTE DE NUEVO Y CAMBIE EL PASO')
     break;
 else
   p=4;
     for(i=1:c)
       fi=feval(f,x(i))
       if(i==1 | i==c)
           suma=suma+fi;
      else
          if(i==p)
              sumav=sumav+(fi*2);
           p=p+3;
       else
         sumav=sumav+(fi*3);
       end 
   end
   end
 end
 suma=suma+sumav; 
   integral=(3*h/8)*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)/80)*(h^4)*f4;
   integral=integral+error