Построение графика функции в турбо-Паскаль

uses graph;
function F(t:real):real;
begin
F:=2*(1+0.7*sin((10e+7)*t)*sin((10e+8)*t)) {функция}
end;
 
var a,b,mx,my,max:real;
    x:real;
    x0,y0,i:integer;
    s:string;
begin
x0:=0;
initgraph(x0,y0,'');{инициализация графического режима}
a:=-1;
b:=1;{интервал по Х}
x0:=getmaxX div 2;{начало координат по Х}
mx:=x0-10;{масштаб по Х}
x:=a; {найдем макс для масштаба по Y}
max:=F(x);
while x<=b do
 begin
  if F(x)>max then max:=F(x);
  x:=x+0.01;
 end;
y0:=getmaxY-50;
my:=(y0-80)/max;
line(0,y0,getmaxX,y0);{оси}
line(x0,60,x0,getmaxY);
for i:=1 to 10 do{максимальное количество засечек в одну сторону}
 begin
  if i<8 then
   begin
    line(x0-3,y0-round(i*my/2),x0+3,y0-round(i*my/2));{засечки на оси У}
    {подпись оси У}
    settextstyle(0,0,1);{горизонтальный шрифт}
    str(i/2:0:1,s);
    outtextXY(x0-35,y0-round(i*my/2),s);{соответственно засечкам}
   end;
  line(x0+round(i*mx/10),y0-3,x0+round(i*mx/10),y0+3); {засечки на оси Х}
  line(x0-round(i*mx/10),y0-3,x0-round(i*mx/10),y0+3);
  {подпись оси Х}
  settextstyle(0,1,1);{вертикальный шрифт}
  str(i/10:0:1,s);
  outtextXY(x0+round(i*mx/10)+5,y0+10,s);
  outtextXY(x0-round(i*mx/10)+5,y0+10,'-'+s);
 end;
settextstyle(0,0,1);{горизонтальный шрифт}
{центр}
outtextXY(x0-15,y0+10,'0');
{подписи концов осей}
outtextXY(getmaxX-20,y0-20,'X');
outtextXY(x0+10,60, 'Y');
{график}
x:=a;
setcolor(12);
moveto(x0+round(x*mx),y0-round(F(x)*my));
while x<=b do
 begin
  lineto(x0+round(x*mx),y0-round(F(x)*my));
  x:=x+0.01;
 end;
{название}
settextstyle(0,0,1);
setcolor(12);
outtextXY(150,10,'y=2*(1+0.7*sin((10e+7)*t)*sin((10e+8)*t))');
outtextXY(250,30,'int.[-1;1]');
readln
end.

uses graph;
const
w1:real=1e+8; {sobstvenni kolebaniy}
tm:real=10e-9; {vremya modelirovaniya}
dt:real=1e-9; {shag modelirovaniya}
n:integer=2; {kol-vo uravneniy}
dk:real=0.0001; {masshtabnyi koef-t}
ti:real=1e-8;{metka}
u1:real=1e+7;
T1:real=0.5e-8;
T2:real=0.2e-8;
type
vec=array[1..2] of real;
matr=array[1..4,1..2] of real;
proc1=procedure(n:integer;var y,dy:vec);
proc2=procedure(n:integer;dt:real;var y,dy:vec;var ier:boolean);
var
i,j,driver,mode,x,z:integer;
ier:boolean;
y,dy:vec;
w2,d,t,f:real;
procedure fct(n:integer;var y,dy:vec) ;far;
begin
f:=sqr(3*(1+0.7*sin(u1*t))*(sin(w1*t)));
dy[1]:=y[2];
dy[2]:=(f-y[1]-y[2]*T2)/(T1*T1);
end;
procedure out(n:integer;dt:real;var y,dy:vec;var ier:boolean); far;
begin
t:=t+dt; {shag vremeny}
if t>tm then ier:=true;
x:=10+round(620/tm*t);
z:=240-round(dk*y[1]);
lineto(x,z);
end;
procedure rk4(n:integer;ier:boolean;h:real;fct:proc1;out:proc2;var y,dy:vec);
var
c:vec;
k:matr;
i,j:integer;
begin
repeat
c:=y;
for i:=1 to 4 do begin
fct (n,y,dy);
for j:=1 to n do begin
k[i,j]:=dy[j]*h;
case i of
3: y[j]:=c[j]+k[i,j];
4: y[j]:=c[j]+(k[1,j]+2*k[2,j]+2*k[3,j]+k[4,j])/6
else y[j]:=c[j]+0.5*k[i,j];
end;end;end;
out (n,h,y,dy,ier);
until ier;
end;
begin
ier:=false;
t:=0; y[1]:=0; y[2]:=0; f:=1/dt;
driver:=detect;
initgraph (driver,mode,'');
line(10,10,10,480);
line(10,240,640,240);
moveto(10,240);
for i:=0 to 10 do begin
x:=10+i*round(620/tm*ti);
line(x,(240-10),x,(240+10));
for j:=1 to 9 do begin
line(x+round(620/tm*ti*0.1*j),(240-5),x+round(620/tm*0.1*ti*j),(240+5));
end;end;
rk4(n,ier,dt,fct,out,y,dy);
readln;
end.

t,f1,f2:array[1..100] of real;