##function derivative (t,y)
##     return exp(t)*y^2+exp(3);
##end

format short

close all

%f = @(x) x^4+2*x^3;

function f = funcao(x)
     f = exp(x);
end

a = 1;
b = 5;

n = 10;
h = (b-a)/n;

method = "simp1/3" % trapz or simp1/3

Int_previous = 10;
Int = 0;

max_Error = 0.01;

%fprintf("k: %0.0f -> t_k: %0.5f & y_k: %0.5f\n",k,t_k,y_k);

%plot(t_k,y_k,'*k');
%hold on;

if strcmp(method,"trapz")

     tic
     do
          for k = 0:n
               if k == 0 || k == n
                    Int = Int + funcao(a+k*h);
               else
                    Int = Int + 2*funcao(a+k*h);
               end
          end

          Int = Int*h/2;

          error = abs(Int-Int_previous);

          Int_previous = Int;

          h = (b-a)/(n++);

     until (error < max_Error)

     disp([Int error n])
     toc

elseif strcmp(method,"simp1/3")
     tic
     do
          for k = 0:n
               if k == 0 || k == n
                    Int = Int + funcao(a+k*h);
               else
                    if rem(k,2) == 0
                         Int = Int + 2*funcao(a+k*h);
                    else
                         Int = Int + 4*funcao(a+k*h);
                    end
               end
          end

          Int = Int*h/3;

          error = abs(Int-Int_previous);

          Int_previous = Int;

          %% SIMPSON ONLY WORKS WITH AN EVEN NUMBER OF INTERVALS
          n++;
          if rem(n,2) ~= 0
               n++;
          endif

          h = (b-a)/n;

     until (error < max_Error)

     disp([Int error n])
     toc

else
     disp('Invalid method!')
     return;
end