Метод Ромберга

clear, clc
 
% подинтегральная ф-ция
f = inline('exp(-x).*sin(5*x).^2');
% пределы
a = 0;
b = 5;
% порядок интегрирования
n = 5;
 
% пустая матрица
R = zeros([n + 1, n + 1]);
% Начальное приближение.  Один интервал по методу трапеций.
R(1,1) = (b - a) / 2 * (f(a) + f(b));
 
% Первый столбец таблицы. 
% Значения интегралов
for i = 1 : n
   h = (b - a) / 2^i;
   s = 0;
   for k = 1 : 2^(i-1)
      s = s + f( a+(2*k-1)*h ) ;
   end
   R(i+1,1) = R(i,1)/2 + h*s;
end
 
% Рассчет уточнений
for j = 1:n
   fac = 1/(4^j-1);
   for m = j:n
      R(m+1, j+1) = R(m+1, j) + fac*(R(m+1,j) - R(m,j));
   end
end
 
R % таблица 
r = R(n+1,1) % итоговый результат
quad(f,a,b) % точное значение для сравнения

R =
    0.0003         0         0         0         0         0
    0.0011    0.0013         0         0         0         0
    0.0012    0.0013    0.0013         0         0         0
    0.0012    0.0012    0.0012    0.0012         0         0
    0.4886    0.6511    0.6944    0.7054    0.7082         0
    0.4916    0.4926    0.4821    0.4787    0.4778    0.4776
 
r =
    0.4916
 
ans =
    0.4918