Численное интегрирование(метод Симпсона)

[(=Hello=)]

Здравствуйте, проверьте программу.
Вычислить значение определенного интеграла с заданной точностью ε.
функция: f(x)=sqrt(1+x^5).
Метод интегрирования Симпсона.

 program test;
uses crt;
var a,b,r,r2,eps:real;
    n:integer;
function f(x:real):real;
begin
f:=sqrt(1+sqr(x)*sqr(x)*x);
end;

function Simpson(a,b:real; n:integer):real;
var s,h:real;
    m,mn:integer;
begin
h:=(b-a)/(2*n);
s:=f(a)+f(b);
mn:=4;
for m:=1 to 2*n-1 do
   begin
   s:=s+mn*f(a+h*m);
   if mn=4 then mn:=2
           else mn:=4;
   end;
   Simpson:=s*h/3;
   end;

begin
writeln ('введите нижнюю границу');
Readln (a);
writeln('введите верхнюю границу');
readln(b);
writeln('введите погрешность');
readln(eps);
n:=2;
r:=simpson(a,b,n);
repeat
 r2:=r;
 n:=n*2;
 r:=simpson(a,b,n);
until abs(r-r2)<eps;
 Writeln;
 Writeln ('Значение интеграла = '  ,r:6:6);
 writeln ('количество разбиений = ',n);
 readkey;
end.

[rubius2008]

Сами проверьте себя:
1. вычислите интеграл "ручками", получите аналитическое решение.
2. потом посмотрите что программа выдает, уменьшите эпсилон, снова посчитайте программой, значение должно получиться ближе к аналитическому. Если не так, то скорей всего что-то в программе не верно.