Вычислить приближенные значения интегралов функций f(x) на интервале [a;b] с точностью 0,001 методом Симпсона - Free Pascal

Формулировка задачи:

Ребятки, помогите пожалуйста с курсачом, задание во вложении

Решение задачи: «Вычислить приближенные значения интегралов функций f(x) на интервале [a;b] с точностью 0,001 методом Симпсона»

Листинг программы

uses Graph,wincrt;
function f(x:real):real;
begin
f:=x*arctan(x)/sqrt(1+x*x);
end;
var a, b, x, y, mx, my, dx: real;
    cX, cY, i, px, py: integer;
    s: string;
    gd, gm : integer;
begin
gd := detect;
initgraph(gd, gm, '');
cX:=getmaxx div 2;
cY:=getmaxy-50;
a:=-5;
b:=5;
mx:=(cX-30)/b;
my:=(cY-50)/f(b);
line(0, cY, getmaxx, cY);
line(getmaxX,cY,getmaxX-10,cY-5);
line(getmaxX,cY,getmaxX-10,cY+5);
line(cX, 0, cX, getmaxy);
line(cX, 0, cX-5, 10);
line(cX, 0, cX+5, 10);
for i:=1 to 15 do
 begin
  line(cX+round(i*mx),cY-3,cX+round(i*mx),cY+3);
  line(cX-round(i*mx),cY-3,cX-round(i*mx),cY+3);
  str(i,s);
  outtextxy(cX+round(i*mx)-3,cY+10,s);
  outtextxy(cX-round(i*mx)-5,cY+10,'-'+s);
  line(cX-3,cY-round(i*my/10),cX+3,cY-round(i*my/10));
  str(i/10:0:1,s);
  outtextxy(cX-30,cY-round(i*my/10),s);
 end;
outtextxy(cX-12,cY+10,'0');
outtextxy(getmaxx-10,cY-20,'X');
outtextxy(cX+90,25, 'Y=x*arctg(x)/sqrt(1+x^2)');
outtextxy(cX+15,5, 'Y');
x:=a;
dx:=0.001;
while x<=b do
 begin
  x:=x+dx;
  px:=cX+round(x*mx);
  y:=f(x);
  py:=cY-round(y*my);
  putpixel(px,py, 8);
 end;
readkey;
end.

Объяснение кода листинга программы

Объявлены переменные a, b, x, y, mx, my, dx, cX, cY, i, px, py, s, gd, gm.
Инициализируются графические режимы и координаты для отрисовки графических примитивов.
Задаются значения переменных a и b для интервала интегрирования.
Вычисляются коэффициенты mx и my для масштабирования осей x и y соответственно.
Строятся оси координат и деления на интервале [a;b] с шагом 1.
В цикле вычисляются значения функции f(x) для каждого деления на интервале [a;b] и ставится точка на соответствующей координате (px,py).
После цикла выводится значение 0 в точку с координатами (cX-12,cY+10).
Строятся оси координат и подписываются их значения.
Вычисляется приближенное значение интеграла функции f(x) на интервале [a;b] с точностью 0,001 и записывается в переменную dx.
В цикле с шагом dx перебираются все значения x на интервале [a;b] и для каждого значения вычисляется соответствующее значение y.
Для каждого значения x и y ставится точка на графике в соответствующей точке с координатами (px,py).
После цикла выводится ключ, который позволяет прервать выполнение программы.