Program diffuse2;
        { define variables...}
const  Ncel = 10;
var   i, j, Timesteps                      : integer;
      Totx, Tott, Delx, Delt, De, mixf     : real;
      Ct1, Ct2, y                          : array[0..Ncel] of real;
      A		                           : array[0..Ncel, 1..3] of real;

begin
        { set variables and initial concentrations...}
  Timesteps := 1;	{ also for Timesteps := 5, 15 }
  De := 1e-5 {cm2/s}; Totx := 1000.0 {cm}; Tott := 100 * 365 * 24 * 60 * 60.0 {s}; Delx := Totx / Ncel;
  for i := 1 to Ncel do Ct2[i] := 0.0; Ct2[0] := 0.56;
	{ calculate the mixing factor ...}
  Delt := Tott / Timesteps;
  mixf := De * Delt / (Delx * Delx);
	{ fill coefficient matrix A ...}
	{ boundary cells ...}
  A[0, 2] := 1; A[0, 3] := 0;
  A[1, 1] := -2 * mixf; A[1, 2] := 1 + 3 * mixf; A[1, 3] := -mixf;
  A[Ncel, 1] := -mixf; A[Ncel, 2] := 1 + mixf;
	{ inner cells ... }
  for i := 2 to Ncel - 1 do
    begin
      A[i, 1] := -mixf; A[i, 2] := 1 + 2 * mixf; A[i, 3] := -mixf;
    end;
	{ decompose A in LU: store L in A[.., 1-2] and U in A[.., 3] ...}
  for i := 1 to Ncel do
    begin
      A[i-1, 3] := A[i-1, 3] / A[i-1, 2];
      A[i, 2] := A[i, 2] - A[i, 1] * A[i-1, 3];
    end;

	{ now diffuse...}
  for i:= 1 to Timesteps do
    begin
      for j := 0 to Ncel do Ct1[j] := Ct2[j];
	{ solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1 ...
	  First, find y in L.y = Ct1 }
      y[0] := Ct1[0];
      for j := 1 to Ncel do y[j] := (Ct1[j] - A[j, 1] * y[j-1]) / A[j, 2];
	{ Now obtain Ct2 in U.Ct2 = y ...}
      Ct2[Ncel] := y[Ncel];
      for j := Ncel - 1 downto 1 do Ct2[j] := y[j] - A[j, 3] * Ct2[j+1];
    end;

  writeln(' Depth below seabottom (m); Cl (mol/L) after', Tott/(31536e3):8:2,' years');
  writeln('             0     ;',Ct2[0]:8:4, Timesteps:5, ' Timesteps');
  for i:=1 to Ncel do
    begin
      writeln('          ', Delx * (i - 0.5) / 100:9:3, Ct2[i]:8:4);
    end;
end.