% rombergtest.m

f = @(x) 4 / (1+x^2); % should converge to pi
r = romberg(f, 0, 1, 5)

% Note: to use a function f1 defined in an m-file f1.m
% use romberg(@f1,0,1,5)

% Try the error function

f = @(x) exp(-x^2);
r = romberg(f, 0, 1, 5)

% Exact value

exact_value = sqrt(pi)*erf(1.0)/2.0