Metoda Romberga - obliczanie całek, implementacja kodu

Witam, potrzebuję pilnie użyć funkcji obliczania całki metodą Romberga w Octave lub Matlab. Znalazłem w sieci kilka takowych funkcji, jednakże nie mam o tym zielonego pojęcia, a fukcje pewnie trzeba jakoś wywołać w głównym programie. Byłby bardzo wdzięczny za (nie oszukujmy się) gotowy kod.
Oto funkcje znalezione przeze mnie:

 function romberg(f,a,b,n)
% Programista Janusz Chojnacki
%Funkcja romberg oblicza całki oznaczone metodą Romberga
fprintf('\n')
disp(' Tablica Romberga')
disp('_________________________________________________________')
disp(' i h Ri,1 Ri,2 Ri,3 ... ')
disp('_________________________________________________________')
h=b-a;
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
fprintf('%2.0f %8.4f %12.4f\n',1,h,R(1,1));
m=1;
for i=1:n-1
 h=h/2;
 S=0;
 for j=1:m
 x=a+h*(2*j-1);
 S=S+feval(f,x);
 end 
4
 R(i+1,1)=R(i,1)/2+h*S;
 fprintf('%2.0f %8.4f %12.4f',i+1,h,R(i+1,1));
 m=2*m;
 for k=1:i
 R(i+1,k+1)=R(i+1,k)+(R(i+1,k)-R(i,k))/(4^k-1);
 fprintf('%12.4f',R(i+1,k+1));
 end
 fprintf('\n');
end 

function Rnum = RombergInt(f,a,b,m)
% Integration algorithm based on Romberg extrapolation
% f - string input for function y = f(x) (e.g. f = 'x.^6')
% a - lower limit of integration
% b - upper limit of integration
% m - maximal number of Romberg iterations
% Rnum - row-vector of numerical approximations for the integral of f(x) between x = a and x = b:
% the entry with index k in Rnum corresponds to the approximation of order O(h^(2k))

R = ones(m,m); % the matrix for Romberg approximations 
              
hmin = (b-a)/2^(m-1); % the minimal step size 

for k = 1 : m
    
    h = 2^(k-1)*hmin; % the step size for the k-th row of the Romberg matrix
    
	x = a : h : b;
	f1 = eval(f);
	k1 = length(f1);
    R(k,1) = 0.5*h*(f1(1) + 2*sum(f1(2:k1-1)) + f1(k1)); % trapezoidal rule for the first column of the Romberg matrix
end
                   
for k = 2 : m % compute the k-th column of the Romberg matrix
    
    for kk = 1 : (m-k+1) % the matrix of Romberg approximations is triangular!
                        
        R(kk,k) = R(kk,k-1)+(R(kk,k-1) - R(kk+1,k-1))/(4^(k-1)-1); % see the Romberg integration algorithm
    end 
        
end

% define the vector Rnum for numerical approximations
Rnum = R(1,:);  

function R = romberg(f,a,b,q)
% ROMBERG Integration using Romberg's method
% Usage: R = ROMBERG(FNAME,A,B,Q)
% R is the Romberg table with Q+2 rows (default Q=2)
% for the integral of F=FNAME(X) from A to B
if nargin<4, q = 2; end 
R = NaN(q+2,q+2);
h = b-a; m = 1;
R(1,1) = h*(feval(f,a)+feval(f,b))/2;
for i = 2:q+2
    h = h/2; m = 2*m;
    R(i,1) = R(i-1,1)/2 + h*sum(feval(f,a+h*(1:2:m)));
    for j = 1:i-1
        R(i,j+1) = R(i,j) + (R(i,j) - R(i-1,j))/(2^(2*j)-1);
    end
end