График функции в turbo pascal

[Trump]

Разработать программу построения графиков по формуле С=sqrt((q*Eo*Es*Nd)/(2*fi-2*U)).
Построить графики для 3 - 5 значений Nd(в диапазоне от 10е18 до 10е23) при изменении U от 0 до –12.В прямоугольной области фиксированного размера построить координатную сетку с оцифровкой.
При изменении значений исходных и расчетных данных масштаб по координатным осям при построении графиков должен изменяться так, чтобы построенные графики максимально полно занимали отведенное для них пространство.
Программу я написал, но работает она не достаточно правильно:
1)графики строятся относительно правильно только, если значения Nd вводить в порядке убывания;
2)даже так я не уверен в правильности графиков;
3)неправильно оцифровывается ось Y.
В каком месте и что нужно исправить?
код программы:

program zk;
uses capwrgr,crt;
var
  Un,Uk,g,n,m,k:integer;
  U,s,Nd:real;
  bul:boolean;
  f:array[1..5] of real;
BEGIN
clrscr;
repeat
writeln('vvedite kol-vo grafikov ot 3 do 5');
readln(k);
if (k=3) or (k=4) or (k=5) then bul:=true else bul:=false;
if bul=false then
begin
clrscr;
writeln('Nevernoe znachenie! Povtorite vvod')end
else writeln;
until bul=true;
for m:=1 to k do
repeat
writeln('vvedite znachenie kontsentratsii Nd',m,' iz diapazona ot 1*10e18 do 1*10e23 m-3');
read(Nd);
if (Nd>=1*10e18) and (Nd<=1*10e23) then bul:=true else bul:=false;
if bul=false then
begin
clrscr;
writeln('Nevernoe znachenie! Povtorite vvod')end
else writeln;
f[m]:=Nd;
until bul=true;
U:=0;
Un:=0;
Uk:=-12;
writeln('vvedite kolichestvo znacheniy napryazheniya');
readln(n);
s:=(Uk-Un)/(n-1);
for m:=1 to k do
  begin
  U:=0;
  for g:=1 to n do
    begin
    a[g]:=U;
    b[m,g]:=cap(f[m],U);
    U:=U+s;
    end;
  end;
wr(n,k);
gr(n,k);
END.

код модуля:

unit capwrgr;
interface
uses crt,graph;
var
  k,p,r,m:integer;
  U:real;
  t:string;
  a:array[0..100] of real;
  b:array[0..100,1..5] of real;
  function cap(var f:array of real;U:real):real;
  procedure wr(var n,k:integer);
  procedure gr(var n,k:integer);
implementation
  function cap(var f:array of real;U:real):real;
  const
  q=1.6e-19;      {elementarny zaryad, Kl}
  Es=11.8;        {otnositel'naya dielektricheskaya pronitsaemost' kremniya}
  Eo=8.85e-12;    {dielektricheskaya postoyannaya}
  fi=0.6;         {poverkhnostny patentsial dlya perehoda zoloto-kremniy}
  begin
    cap:=sqrt((q*Eo*Es*f[m])/(2*(fi-U)));
  end;
  procedure wr;
  var
    j,m:integer;
  begin
  clrscr;
    for m:=1 to k do
      begin
        writeln('Napryazhenie U,V           |    Yomkost C',m,',F/m2');
        for j:=1 to n do
          writeln(a[j],'    |',b[m,j]);
          readln;
      end;
  readkey;
  end;
  procedure gr;
   var
   i,j,
   x,y,
   grdriver,
   grmode:integer;
begin
  grdriver:=vga;
  grmode:=vgahi;
  initgraph(grdriver, grmode, 'd:\bp\bgi');
  for m:=1 to k do
    begin
    setlinestyle (userbitln, $8888,normwidth);
    for i:=0 to n do line (getmaxx-round(580/(n+1)),450-round((420/(n-1))*i),
    getmaxx-round((580*(n+1))/(n+1)),450-round((420/(n-1))*i));{gor}
    setlinestyle (userbitln, $8888, normwidth);
    for i:=0 to n do line (60+round((580*i)/n),30,
    60+round((580*i)/n),450);{ver}
    for i:=0 to n do
    begin
    str(round(a[n-i]),t);
    outtextxy(55+round((580*i)/n+1),455,t);end;{otsifrovka X}
    for i:=0 to n do
    begin
   { str(((b[m,n-i])*1e7):7:2,t);}
    str((((i+1)*(b[m,n-n]-b[m,n+1])/(n+1))*1e7):7:2,t);
    outtextxy(0,446-round((420/(n-1))*i),t);end;{otsifrovka Y}
    setlinestyle(0,1,1);
    moveto(getmaxx-round(580/(n+1)),30);
    for j:=1 to n do
      begin
        x:=getmaxx-round((580*j)/(n));
        y:=450-round((420*(b[m,j]-b[m,n]))/(b[m,1]-b[m,n]));
        setcolor(m);
        lineto(x,y);
        circle(x,y,3);
      end;
    setcolor(15);
    outtextxy(0,470,'dlya vyhoda nazhmite "enter"');
    outtextxy(0,0,'C*1e7, f');
    outtextxy(600,445,'U, V');
    readln;
    cleardevice;
    end;
  closegraph;
end;
end.