## Author: M. Turhan ÇOBAN
## Created: 2021-06-13

function retval = newton_cotes_integral (f,a,b,n,n_formula)
  % f f(x) function
  % a lower integration limit
  % b upper integration limit
  % n number of division ao integration area
  % n_formula:
  % 2 : Trapezoidal rule
  % 3 : Simpson 1/3 rule
  % 4 : Simpson 3/8 rule
  % 5 : Boole rule
  retval=0.0;
if n_formula==2
  h=(b-a)/n;
  ff=f(a)+f(b);
  for i=1:n-1
    x=a+i*h;
    ff=ff+2.0*f(x);
   endfor 
   ff=h/2.0*ff;
   retval =ff;
elseif n_formula==3
  h=(b-a)/n;
  h1=h/2.0;
  ff=0.0;
  fi=zeros(1,n_formula);
  for p=1:n
    for i=1:n_formula
      x=a+(i-1)*h1+(p-1)*h;
      fi(i)=f(x);
    endfor
   ff=ff+h*(fi(1)+4.0*fi(2)+fi(3))/6.0;
   endfor 
   retval =ff;
elseif n_formula==4
  h=(b-a)/n;
  h1=h/(n_formula-1.0);
  ff=0.0;
  fi=zeros(1,n_formula);
  for p=1:n
    for i=1:n_formula
      x=a+(i-1)*h1+(p-1)*h;
      fi(i)=f(x);
    endfor
   ff=ff+h*(fi(1)+3.0*fi(2)+3.0*fi(3)+fi(4))/8.0;
   endfor 
   retval =ff;
elseif n_formula==5
  h=(b-a)/n;
  h1=h/(n_formula-1.0);
  ff=0.0;
  fi=zeros(1,n_formula);
  for p=1:n
    for i=1:n_formula
      x=a+(i-1)*h1+(p-1)*h;
      fi(i)=f(x);
    endfor
   ff=ff+h*(7.0*fi(1)+32.0*fi(2)+12.0*fi(3)+32.0*fi(4)+7.0*fi(5))/90.0;
   endfor 
   retval =ff; 
endif
endfunction
