[hatabi4]

Добрый вечер всем.
Прошу помощи у знающих людей. По предмету вычислительная математика дали задание, которое никак не получается. Конец сессии нужно срочно сдать. суть вот в чем:
Нужно провести интерполяцию функции F=x*x+4*x+3, используя интерполяционную формулу Ньютона(метод разделенных разностей) на интервале [-10..10] взяв количество узлов интерполяции равных 3. Построить графики исходной функции и интерполирующей функции, а так же график погрешности интерполяции.
Сделать нужно, используя вот эти формулы:


В моем случае уравнение получается вида:
Pn(x)=f(x[0])+f(x[i],x[i+1])*(x-x[0])+f(x[i],x[i+1],x[i+2])*(x-x[0])*(x-x[1]).

Написал функцию, но работать она не хочет:

Код:

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, TeEngine, Series, ExtCtrls, TeeProcs, Chart, Math;

type
  TForm1 = class(TForm)
    Chart1: TChart;
    Chart2: TChart;
    Series1: TLineSeries;
    Series2: TLineSeries;
    Series3: TLineSeries;
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

const
  xmin=-10;
  xmax=10;
  n=3;

var
  Form1: TForm1;
  i,j: integer;
  x1,x2,x3,h: real;
  x,fx: array [1..3] of real;

implementation

{$R *.dfm}

function yf(x:real):real;
begin
  yf:=x*x+4*x+3;
end;


Function Newton(n: integer; d:real; x,y :array of real):real;
var
  l:real;
  j,i,k:integer;
  p,p2: real;
  y1,y2: array of real;
begin
L:=fx[1];
P:=1;
p2:=1;
for j:=1 to n do
begin
  P:=P*(D-X[j-1]);
    for i:=0 to (n-j) do // <- ПЕРВАЯ РАЗНОСТЬ
    begin
      Y[i]:=(y[i+1]-y[i])/(x[i+j]-x[i]);  
    end;
  p2:=p2*(d-x[j-1])*(d-x[j]); 
    for k := 1 to (n-j) do    // <- ВТОРАЯ РАЗНОСТЬ, ЗДЕСЬ ВСЯ "СОЛЬ" 
      begin
        y[k]:=(((y[k+2]-y[k+1])/(x[k+j]-x[k]))-y1[k])/(x[k+2]-x[k]);  
      end;
  L:=L+P*y[i]+p2*y[k];
end;
  Newton:=-l;
end;



procedure TForm1.FormCreate(Sender: TObject);
begin

  x1:=xmin;
  while x1<=xmax do
  begin
    Series1.AddXY(x1,yf(x1));
    x1:=x1+1;
  end;

  x2:=-1;
  for i:=1 to n do
  begin
    x[i]:=x2;
    fx[i]:=yf(x[i]);
    x2:=x2+1;
  end;


  x3:=xmin;
  while x3<=xmax do
  begin
    for i:=1 to 10 do
    begin
    Series2.AddXY(x3,Newton(n,x3,x,fx));
    Series3.AddXY(x3, abs (yf(x3)-Newton(n,x3,x,fx)) );
    x3:=x3+0.01;
    end;
  end;

end;

end.

Математический я все понимаю как делать и на листике оно то получается, а программно никак. Помогите пожалуйста с кодом функции.