library IntervalMilneMethod_ExamplaryLibrary;
uses System.SysUtils, System.Classes, IntervalArithmetic32and64;
{$R *.res}

// The library includes all functions needed for an interval version of Milne's
// multistep method determined in IntervalMultistepMethods program.
// The functions defined in this library concerns the initial value problem
// of the form
//                  y' = (y - t)/(y + t), y(0) = 4,
// Note: Do not use this unit for any Delphi's compiler older than XE!

function FTY (const T, Y : interval) : interval;
begin
  Result:=(Y-T)/(Y+T)
end {FTY};

function D1FDT (const T, Y : interval) : interval;
var YT2 : interval;
begin
  YT2:=Y+T;
  YT2:=YT2*YT2;
  Result:=-2*Y/YT2
end {D1FDT};

function D1FDY (const T, Y : interval) : interval;
var YT2 : interval;
begin
  YT2:=Y+T;
  YT2:=YT2*YT2;
  Result:=2*T/YT2
end {D1FDY};

function D2FDT2 (const T, Y : interval) : interval;
var YT3 : interval;
begin
  YT3:=Y+T;
  YT3:=YT3*YT3*YT3;
  Result:=4*Y/YT3
end {D2FDT2};

function D2FDTDY (const T, Y : interval) : interval;
var YT3 : interval;
begin
  YT3:=Y+T;
  YT3:=YT3*YT3*YT3;
  Result:=2*(Y-T)/YT3
end {D2FDTDY};

function D2FDY2 (const T, Y : interval) : interval;
var YT3 : interval;
begin
  YT3:=Y+T;
  YT3:=YT3*YT3*YT3;
  Result:=-4*T/YT3
end {D2FDY2};

function D3FDT3 (const T, Y : interval) : interval;
var YT4 : interval;
begin
  YT4:=Y+T;
  YT4:=YT4*YT4;
  YT4:=YT4*YT4;
  Result:=-12*Y/YT4
end {D3FDT3};

function D3FDT2DY (const T, Y : interval) : interval;
var YT4 : interval;
begin
  YT4:=Y+T;
  YT4:=YT4*YT4;
  YT4:=YT4*YT4;
  Result:=4*(T-2*Y)/YT4
end {D3FDT2DY};

function D3FDTDY2 (const T, Y : interval) : interval;
var YT4 : interval;
begin
  YT4:=Y+T;
  YT4:=YT4*YT4;
  YT4:=YT4*YT4;
  Result:=4*(2*T-Y)/YT4
end {D3FDTDY2};

function D3FDY3 (const T, Y : interval) : interval;
var YT4 : interval;
begin
  YT4:=Y+T;
  YT4:=YT4*YT4;
  YT4:=YT4*YT4;
  Result:=12*T/YT4
end {D3FDY3};

function D4FDT4 (const T, Y : interval) : interval;
var YT5 : interval;
begin
  YT5:=Y+T;
  YT5:=YT5*YT5;
  YT5:=YT5*YT5;
  YT5:=YT5*(Y+T);
  Result:=48*Y/YT5
end {D4FDT4};

function D4FDT3DY (const T, Y : interval) : interval;
var YT5 : interval;
begin
  YT5:=Y+T;
  YT5:=YT5*YT5;
  YT5:=YT5*YT5;
  YT5:=YT5*(Y+T);
  Result:=12*(3*Y-T)/YT5
end {D4FDT3DY};

function D4FDT2DY2 (const T, Y : interval) : interval;
var YT5 : interval;
begin
  YT5:=Y+T;
  YT5:=YT5*YT5;
  YT5:=YT5*YT5;
  YT5:=YT5*(Y+T);
  Result:=24*(Y-T)/YT5
end {D4FDT2DY2};

function D4FDTDY3 (const T, Y : interval) : interval;
var YT5 : interval;
begin
  YT5:=Y+T;
  YT5:=YT5*YT5;
  YT5:=YT5*YT5;
  YT5:=YT5*(Y+T);
  Result:=12*(Y-3*T)/YT5
end {D4FDTDY3};

function D4FDY4 (const T, Y : interval) : interval;
var YT5 : interval;
begin
  YT5:=Y+T;
  YT5:=YT5*YT5;
  YT5:=YT5*YT5;
  YT5:=YT5*(Y+T);
  Result:=-48*T/YT5
end {D4FDY4};

function D1FTY (const T, Y : interval) : interval;
begin
  Result:=D1FDT(T,Y)+D1FDY(T,Y)*FTY(T,Y)
end {D1FTY};

function D2FTY (const T, Y : interval) : interval;
var F, F2 : interval;
begin
  F:=FTY(T,Y);
  F2:=(D2FDY2(T,Y)*F+2*D2FDTDY(T,Y))*F;
  Result:=D2FDT2(T,Y)+F2+D1FDY(T,Y)*D1FTY(T,Y)
end {D2FTY};

function D3FTY (const T, Y : interval) : interval;
var F, F2, F3 : interval;
begin
  F:=FTY(T,Y);
  F3:=(D3FDY3(T,Y)*F+3*D3FDTDY2(T,Y))*F;
  F3:=(F3+3*D3FDT2DY(T,Y))*F;
  F3:=D3FDT3(T,Y)+F3;
  F2:=(D2FDTDY(T,Y)+D2FDY2(T,Y)*F)*D1FTY(T,Y);
  Result:=F3+3*F2+D1FDY(T,Y)*D2FTY(T,Y);
end {D3FTY};

function PSITY (const T, Y : interval) : interval;
var DF, F, F2, F3, F4 : interval;
begin
  F:=FTY(T,Y);
  F4:=(D4FDY4(T,Y)*F+4*D4FDTDY3(T,Y))*F;
  F4:=(F4+6*D4FDT2DY2(T,Y))*F;
  F4:=(F4+4*D4FDT3DY(T,Y))*F;
  F4:=D4FDT4(T,Y)+F4;
  F3:=(D3FDY3(T,Y)*F+2*D3FDTDY2(T,Y))*F;
  F3:=6*(D3FDT2DY(T,Y)+F3);
  DF:=D1FTY(T,Y);
  F2:=D2FDY2(T,Y);
  F3:=(F3+3*F2*DF)*DF;
  F2:=4*(D2FDTDY(T,Y)+F2*F)*D2FTY(T,Y);
  F3:=F3+F2+D1FDY(T,Y)*D3FTY(T,Y);
  Result:=F4+F3
end {PSI};

exports FTY, PSITY;

begin
end.