close all;
clear

 fun = @(x) x.^(0.3)
   
  
% definition av integration intervalet [n,p]

n= 0;
p= 1.0

% gammalt  matlab-funktion
 Q1 =  quadl(fun,n,p);

Q =  integral(fun,n,p);

fun_exact = 0;
x_exact = linspace(n, p);
N = size(x_exact,2)

  fun_calc = 0;
m= (n+p)/2.0;
x_calc = [n,m,p];
N_calc = 3;

% Simpsons formel

int_simpson  = ((p-n)/6.0)*(fun(n) +  4.0*fun(m) +  fun(p));
  
  for i = 1:N
	    func_exact(i) =fun(x_exact(i));
end


for i = 1:N_calc
	  frac1 = ((x_calc(i)- m)*(x_calc(i)- p))/((n - m)*(n-p));

          frac2 = ((x_calc(i)- n)*(x_calc(i)- p))/(( m - n)*(m  - p));
	    
          frac3 = ((x_calc(i)- n)*(x_calc(i)- m))/((p-n)*(p - m));
	    
         fun_calc(i) = fun(n)*frac1 + fun(m)*frac2 + fun(p)*frac3;

end

% trapetsmetoden
int_trap = trapz(x_calc,  fun_calc);


figure
plot(x_exact,func_exact, ' r-',  'LineWidth',2)

  
xverts = [x_calc(1:end-1);x_calc(1:end-1);x_calc(2:end);x_calc(2:end)];
yverts = [zeros(1,N_calc-1);fun_calc(1:end-1);fun_calc(2:end);zeros(1,N_calc-1)];
hold on


patch(xverts, yverts,'g','LineWidth',1.5)

plot(x_calc, fun_calc ,'b','LineWidth',3);

error1 = abs(Q-int_simpson)
error2 = abs(Q-int_trap)
  
  xlabel('x')
  ylabel('funktion f(x)= x^{0.3}')
  legend('exact ','Area (integral)','approximation ');

           title(['Trapetsmetoden:',num2str(int_trap), 'Adaptive integral:',num2str(Q),'error:',num2str(error2)]);