Решение явной четырехточечной схемой методом конечных разностей

//uses crt;
const
  mf = 500;
 
type
  vector = array[0..mf] of real;
 
var{раздел описания переменных, которые мы будем использовать в
  программе}
  i, N: integer;
  T, TT, ro, roro, u, uu, P, PP: vector;
  h, tau, beta, lambda: real;
  T1, T0, Tr, ro1, ro0, ror, u1, u0, ur, P1, Pr, L, t_end, time, c, a, q, P0, k, Tz, Ty: real;
  f, g: text;
 
begin
  //clrscr;
  
  L := 0.05;
  
  {с клавиатуры вводим все необходимые входные параметры}
  Writeln('Введите количество узлов по пространственной координате,N');
  //N := 10;
  Readln(N);
  //Writeln('Введите окончание по времени, t_end');
  //Readln(t_end); 
  t_end := 480;
  
  Writeln('Производим расчет');
  
  
  
  {определяем расчетный шаг сетки по пространственной координате}
  h := L / (N - 1); 
  
  {начальная скорость}
  u1 := 0.0;
  {"нулевая" температура}
  T0 := 293.15;
  //T0:=20;
  {ро первая}
  ro1 := 998.23;
  {тепловой поток}
  q := -10000;
  {коэффициент объемного расширения}
  beta := 1.82 / 10000;
  {теплоемкость}
  c := 4200;
  {плотность нулевая}
  ro0 := 998.23;
  {теплопроводность}
  lambda := 0.599;
  {давление первое}
  P1 := 0.0;
  {скорость первая}
  u1 := 0.0;
  {скорость на границе}
  ur := 0.0;
  {давление на границе}
  Pr := 0.0;
  {плотность на границе}
  ror := 998.23;
  {первое значение температуры после отчета времени и начала шага}
  T1 := (q * h - T0) * -1;
  Tr := T0;
  {значения давления на глубине}
  
  
  {определяем коэффициент температуропроводности}
  a := lambda / (ror * c);
  {определяем расчетный шаг сетки по времени, руководствуясь
  условием устойчивости}
  //tau := ((ro0 * c * sqr(h)) /( 2 * lambda))*0.25;
  tau:=0.001;
  Writeln('цикл 1');
  
  ro[0] := ro1;
  T[0] := T1;
  u[0] := u1;
  P[0] := ro1 * 9.8 * (L - h);
  {определяем поле температуры в начальный момент времени}
  for i := 1 to N - 2 do
  begin
    T[i] := T0;
    ro[i] := ro0;
    u[i] := u0;
    P[i] := ro0 * 9.81 * (L - (h * i));
    //  print(P[i])
  end;
  {определяем значения температуры на границе}
  
  T[N - 1] := Tr;
  ro[N - 1] := ror;
  u[N - 1] := u0;
  P[N - 1] := 0;
  
  {задаем коэффициент k}
  k := 2193361.0;
  
  {проводим интегрирование слау}
  time := 0;
  
  Assign(g, 'E:\kursovaya\pressure.txt'); 
  Rewrite(g);
  writeln(g, 'Distance'#9'temperature'#9'pressure'#9'density'#9'speed'#9'time');
  
  while time < t_end do {используем цикл с предусловием}
  begin
    time := time + tau;
    {запоминаем поле температуры на предыдущем слое по времени}
    Writeln('цикл 2');
    Writeln('Time:' + time);
    Writeln('цикл расчета');
    for i := 0 to N - 1 do
    begin
      if i = 0 then
      begin
        TT[i] := T[i];
        TT[0] := (q * h - T[1]) * -1;
        P[0]:=P[1];
      end
      else
        TT[i] := T[i];
          TT[0] := (q * h - T[1]) * -1;
      uu[i] := u[i];
      roro[i] := ro[i];
    end;   
    
    
    
    for i := 1 to N - 2 do
    begin
     if time>= 2*tau then
      begin
     
      PP[i] := (roro[i] - ro0 * (1.0 - beta * (TT[i] - T0))) * k;
       WriteLn(PP[i]);
       end;
      u[i] := ((-1 / roro[i]) * ((P[i + 1] - P[i]) / h) - uu[i] * ((uu[i + 1] - uu[i]) / h)) * tau + uu[i];
      ro[i] := -(roro[i + 1] * uu[i + 1] - roro[i] * uu[i]) / h * tau + roro[i];
      Tz := (TT[i + 1] - 2.0 * TT[i] + TT[i - 1]) / sqr(h);
      Ty := (TT[i + 1] - TT[i]) / h;
      T[i] := (a * Tz - uu[i] * Ty) * tau + TT[i];
      if i=N-2 then 
      begin
      a:=lambda/ro[i]*c;
      end;
       end;
 
    
    
     
    writeln(g, (h * i):6:3, #9, T[i]:6:3, #9, PP[i]:6:3, #9, ro[i]:6:3, #9, u[i]:6:3, #9, time:6:3);
  end;
  close(g);
end.