%          f_1     f_2          f_3
f  = @(x) [sin(x); 0.01*x.^2;   atan(x)];
a  = -1;
b  =  10;
N  = [2 5 10 20 100];   nN = size(N,2);
nn = 100;
xx = a:(b-a)/nn:b;

I_o = zeros(3,nN+1);   I_l = I_o;   I_s = I_o;

% Obdelnikove    pravidlo
I_o(:,1) = ( b - a )   * (          f((a+b)/2)        );
% Lichobeznikove pravidlo
I_l(:,1) = ( b - a )/2 * ( f(a)                + f(b) );
% Simpsonovo     pravidlo
I_s(:,1) = ( b - a )/6 * ( f(a) + 4*f((a+b)/2) + f(b) );

for i = 1:nN
  for j = 1:N(i)
    a_ = a + (j-1)*(b-a)/N(i);
    b_ = a + (j  )*(b-a)/N(i);
    c_ = (a_ + b_)/2;
    I_o(:,i+1) = I_o(:,i+1) + ( b_ - a_ )   * (           f(c_)       );
    I_l(:,i+1) = I_l(:,i+1) + ( b_ - a_ )/2 * ( f(a_)          + f(b_) );
    I_s(:,i+1) = I_s(:,i+1) + ( b_ - a_ )/6 * ( f(a_) + 4*f(c_)+ f(b_) );
  end
end

disp(I_o);
disp(I_l);
disp(I_s);

% plot
figure(1);   hold on;
yy = f(xx);
plot(xx,yy(1,:),'-r'); % f_1
plot(xx,yy(2,:),'-g'); % f_2
plot(xx,yy(3,:),'-b'); % f_3

