Написать программу вычисления первой и второй производной табличной функции с помощью полинома Ньютона или Лагранжа

{$APPTYPE CONSOLE}
uses sysutils;
const max=50;
type
        TInterval=array[1..2,1..max] of real; {Сетка интерполяции}
    {1 строка - значения узлов в сетке интерполяции
     2 строка - значения функции в узлах интерполяции}
 
        TIntervalOutput=array[1..4,1..max] of real;
    {1 строка - значения узлов в сетке интерполяции
     2 строка - значения функции в узлах интерполяции
     3 строка - значения первой производной функции
     4 строка - значения второй производной функции}
 
 
function f(x:real):real;
{Основная функция}
 begin
 f:=exp(x*x);
 end;
 
 
procedure Newton(x:TInterval; N1,N2:integer; A1,B1:real; var Y:TIntervalOutput);
{Метод Ньютона для интерполирования}
var
    g,h:real;
    i:integer;
 
    function New1(i:integer):real;
    {вычисление очередного члена первой формулы Ньютона - рекурсивна
     вход - порядковый номер}
    begin
        if i=1 then New1:=g else New1:=(New1(i-1)*(g-i+1))/i;
    end;
 
    function d(i,n:integer):real;
    {вычисление конечных разностей - рекурсивна
     вход - номер строки и столбца в таблице разностей}
    begin
        if i=1 then d:=X[2,n+1]-X[2,n] else d:=d(i-1,n+1)-d(i-1,n);
    end;
 
    function newton1(UzelX:real):real;
    {реализация первой формулы Ньютона}
    var
        Next,t1,t2,s:real; {очередной член ф. Ньютона}
        i,n:integer;
    begin
        i:=1;
        n:=trunc((UzelX-x[1,1])/(x[1,2]-x[1,1]))+1; {номер строки таблицы}
        S:=X[2,n];  {сумма членов}
        g:=(UzelX-x[1,n])/(x[1,2]-x[1,1]);
        repeat
            t1:=New1(i);
            t2:=d(i,n);
            Next:=t1*t2;
            S:=S+Next;
            inc(i);
        until (i>(N1-n));  {до конца таблицы}
        Newton1:=S;
    end;
 
    function New2(i:integer):real;
    {вычисление очередного члена второй формулы Ньютона - рекурсивна
     вход - порядковый номер}
    begin
        if i=1 then New2:=g else New2:=(New2(i-1)*(g+i-1))/i;
    end;
 
    function newton2(UzelX:real):real;
    {реализация второй формулы Ньютона}
    var
        Next,t1,t2,s:real;   {очередной член ф. Ньютона}
        n,i:integer;
    begin
        n:=trunc((UzelX-X[1,1])/(x[1,2]-x[1,1])+2);    {номер строки таблицы}
        if n>n1 then n:=n1;  {если за пределами таблицы - экстраполяция}
        g:=(UzelX-X[1,n])/(x[1,2]-x[1,1]);
        S:=x[2,n];  {сумма членов}
        i:=1;
        repeat
            t1:=New2(i);
            t2:=d(i,n-i+1);
            Next:=t1*t2;
            S:=S+Next;
            inc(i);
        until (i>n);  {до конца табл}
        Newton2:=S;
    end;
 
begin
    h:=(B1-A1)/N2;
    for i:=1 to N2 do
    begin
        y[1,i]:=A1+h*(i-1);  {Заполнение таблицы значений}
        if y[1,i]>((x[1,N1]+x[1,1])/2) then y[2,i]:=newton2(y[1,i])
            else y[2,i]:=newton1(y[1,i]);
    end;
 
    {Вычисление первой производной с помощью конечно-разностных отношений}
    y[3,1]:=(y[2,2]-y[2,1])/h;
    for i:=2 to N2-1 do y[3,i]:=(y[2,i+1]-y[2,i-1])/(2*h);
    y[3,N2]:=(y[2,N2]-y[2,N2-1])/h;
 
    {Вычисление второй производной с помощью конечно-разностных отношений}
    for i:=2 to N2-1 do y[4,i]:=(y[2,i+1]-2*y[2,i]+y[2,i-1])/(h*h);
end;
 
 
procedure Lagrange(x:TInterval; N1,N2:integer; A1,B1:real; var Y:TIntervalOutput);
{Интерполяция по формуле Лагранжа}
 
    function L(i,j: integer; UzelX:real):real;
    {Функция вычисления очередного коэффициента формулы Лагранжа по схеме Эйткена - рекурсивна}
    begin
        if j=0 then
    L:=(x[2,i]*(x[1,i+1]-UzelX)-X[2,i+1]*(X[1,i]-UzelX))/(X[1,i+1]-X[1,i])
        else
    L:=(L(i,j-1,UzelX)*(X[1,i+j+1]-UzelX)-L(i+1,j-1,UzelX)*(X[1,i]-UzelX))/(X[1,i+j+1]-X[1,i]);
    end;
 
var
    h:real;
    i:integer;
begin
    h:=(B1-A1)/N2;
    for i:=1 to N2 do
    begin
        y[1,i]:=A1+h*(i-1);
        y[2,i]:=L(1,N1-1,y[1,i])
    end;
 
    {Вычисление первой производной с помощью конечно-разностных отношений}
    y[3,1]:=(y[2,2]-y[2,1])/h;
    for i:=2 to N2-1 do y[3,i]:=(y[2,i+1]-y[2,i-1])/(2*h);
    y[3,N2]:=(y[2,N2]-y[2,N2-1])/h;
 
    {Вычисление второй производной с помощью конечно-разностных отношений}
    for i:=2 to N2-1 do y[4,i]:=(y[2,i+1]-2*y[2,i]+y[2,i-1])/(h*h);
end;
 
var
    A,B,A1,B1,h:real;
    N1,N2,Method,i:integer;
    x:TInterval;
    y:TIntervalOutput;
begin
    Writeln;
    Write('Vvedite nachalo intervala f(x): ');
    Readln(A);
    Write('Vvedite konec intervala f(x): ');
    Readln(B);
    Write('Vvedite kolichestvo uzlov: ');
    Readln(N1);
    Write('Vvedite nachalo intervala interpolyacii: ');
    Readln(A1);
    Write('Vvedite konec intervala interpolyacii: ');
    Readln(B1);
    Write('Vvedite kolichestvo uzlov intervala interpolyacii: ');
    Readln(N2);
    Write('Metod (1-Newton, 2-Lagrange): ');
    Readln(Method);
  if Method=2 then A:=A+0.001;
    h:=(B-A)/N1;
    for i:=1 to N1 do
    begin
        x[1,i]:=A+h*(i-1);
        x[2,i]:=f(x[1,i]);
    end;
    case Method of
        1: Newton(x,N1,N2,A1,B1,Y);
        2: Lagrange(x,N1,N2,A1,B1,Y);
    end;
 
    {вывод результатов}
    Writeln('Metod: ', Method);
    for i:=1 to N2 do
        writeln('x',i:2,'=',y[1,i]:2:4,'  y',i:2,'=',y[2,i]:2:4,'  y',#39,i:2,'=',y[3,i]:2:4,'  y',#39#39,i:2,'=',y[4,i]:2:4);
    writeln;
    for i:=1 to N1 do
        writeln('x',i:2,'=',x[1,i]:2:4,'  y',i:2,'=',x[2,i]:2:4);
    readln;
end.