График изотермы Ван-дер-Ваальса

[Dfds]

Используется синтаксис Delphi

unit Unit1;
 
{$mode objfpc}{$H+}
 
interface
 
uses
  Classes, SysUtils, FileUtil, TAGraph, TASeries, Forms, Controls, Graphics,
  Dialogs, StdCtrls;
 
const
   R=8.314;
   k=1000;
   Zc=0.375;
 
type
 
  { TForm1 }
 
  TForm1 = class(TForm)
    Button1: TButton;
    Button2: TButton;
    Chart1: TChart;
    Chart1LineSeries1: TLineSeries;
    Edit6: TEdit;
    Edit7: TEdit;
    Memo1: TMemo;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
  private
    { private declarations }
  public
    { public declarations }
  end;
 
var
  Form1: TForm1;
  nu,P,T,V,Vkon,Vnach,dV,a,b,Tk,Vk,Pk: real;
  i: integer;
 
implementation
 
{$R *.lfm}
 
{ TForm1 }
 
procedure TForm1.Button1Click(Sender: TObject);
begin
memo1.clear;
Vnach:=strtofloat(edit6.text);
Vkon:=strtofloat(edit7.text);
//воздух
Tk:=132.3;
Pk:=3640000;
//воздух
T:=0.8*Tk;
nu:=1;
Vk:=(nu*Zc*R*Tk)/(Pk);
a:=(27*sqr(Tk)*sqr(R))/(64*Pk);
b:=(Tk*R)/(8*Pk);
memo1.lines.add('b =' + floattostr((b)));
memo1.lines.add('a =' + floattostr((a)));
memo1.lines.add('Vk =' + floattostr((Vk)));
end;
 
procedure TForm1.Button2Click(Sender: TObject);
begin
dV:=(Vkon-Vnach)/(k);
//воздух
Tk:=132.3;
Pk:=3640000;
//воздух
T:=0.8*Tk;
for i:=0 to k do begin
V:=Vnach+dV*i;
P:=(nu*R*T)/(V-nu*b)-(a*nu*nu)/(V*V);
chart1.title.visible:=true;
chart1.title.text.text:='Ван-дер-Ваальс';
chart1lineseries1.addXY(V,P,'',clRed);
end;
end;
end.

Помогите, пожалуйста, найти ошибку
Хочу получить изотерму, как на

Получаю

Иногда выходила нужная изотерма, когда Vnach=1, Vkon=0.0001, но там выходило отрицательное давление