PROGRAM IME41;
CONST
   np =  2;
   mp =  1;
TYPE
   RealArrayNPbyNP = ARRAY [1..np,1..np] OF real;
   RealArrayNPbyMP = ARRAY [1..np,1..mp] OF real;
   IntegerArrayNP = ARRAY [1..np] OF integer;
PROCEDURE gaussj(VAR a: RealArrayNPbyNP;
                     n: integer;
                 VAR b: RealArrayNPbyMP;
                     m: integer);
{
 Resi system linearnich rovnic s mnoha pravymi stranami
 Gauss-Jordanovou eliminaci s uplnym hledanim hlavniho prvku.
 Matice a[n,n], m pravych stran v matici b[n,m]
 na vystupu je v a inverzni matice a^-1
 a v b jsou reseni rovnic
}
VAR
   big,dum,pivinv: real;
   i,icol,irow,j,k,l,ll: integer;
   indxc,indxr,ipiv: ^IntegerArrayNP;
BEGIN
   new(indxc);
   new(indxr);
   new(ipiv);
   FOR j := 1 TO n DO ipiv^[j] := 0;
   FOR i := 1 TO n DO BEGIN
      big := 0.0;
      FOR j := 1 TO n DO
         IF ipiv^[j] <> 1 THEN
            FOR k := 1 TO n DO
               IF ipiv^[k] = 0 THEN
                  IF abs(a[j,k]) >= big THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               ELSE IF ipiv^[k] > 1 THEN BEGIN
                  writeln('pause 1 in GAUSSJ - singular matrix');
                  readln
               END;
      ipiv^[icol] := ipiv^[icol]+1;
      IF irow <> icol THEN BEGIN
         FOR l := 1 TO n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 TO m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      END;
      indxr^[i] := irow;
      indxc^[i] := icol;
      IF a[icol,icol] = 0.0 THEN BEGIN
         writeln('pause 2 in GAUSSJ - singular matrix');
         readln
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 TO n DO
         a[icol,l] := a[icol,l]*pivinv;
      FOR l := 1 TO m DO
         b[icol,l] := b[icol,l]*pivinv;
      FOR ll := 1 TO n DO
         IF ll <> icol THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 TO n DO
               a[ll,l] := a[ll,l]-a[icol,l]*dum;
            FOR l := 1 TO m DO
               b[ll,l] := b[ll,l]-b[icol,l]*dum
         END
   END;
   FOR l := n DOWNTO 1 DO
      IF indxr^[l] <> indxc^[l] THEN
         FOR k := 1 TO n DO BEGIN
            dum := a[k,indxr^[l]];
            a[k,indxr^[l]] := a[k,indxc^[l]];
            a[k,indxc^[l]] := dum
         END;
   dispose(ipiv);
   dispose(indxr);
   dispose(indxc)
END;

   Function mexp(x:real):real;
BEGIN
   if(x > 0) then writeln('x>0',x);
   if(x > -84) then
      mexp := exp(x)
   else
      mexp := 0.;
END;


VAR a1,a2,b1,b2,h,x,xmax,hu,hv:real;
Var n,nkrok,npkrok:integer;
VAR mate:RealArrayNPbyNP;
    vecb:RealArrayNPbyMP;
    yold,ynew,ypom:array [1..2] of real;

BEGIN
    a1:=998; a2:=999; b1:=1998; b2:=1999;
    x:=0; yold[1]:=1;yold[2]:=0;
    writeln('zadej xmax, nkrok, npkrok');
    readln(xmax,nkrok,npkrok);
    h := xmax/nkrok;
    mate[1,1]:= 1-a1*h;
    mate[1,2]:= -b1*h;
    mate[2,1]:= a2*h;
    mate[2,2]:= 1+b2*h;
    vecb[1,1]:=1;vecb[2,1]:=1;
    gaussj(mate,2,vecb,1);
    writeln(x:14,yold[1]:14,yold[2]:14);
    for n:=1 to nkrok do begin
        ynew[1] := mate[1,1]*yold[1]+mate[1,2]*yold[2];
        ynew[2] := mate[2,1]*yold[1]+mate[2,2]*yold[2];
        x:=x+h;
        yold[1]:=ynew[1];
        yold[2]:=ynew[2];
         hu := 2*exp(-x) - mexp(-1000*x);
         hv := -exp(-x) + mexp(-1000*x);
        if (n=nkrok) or (n mod npkrok = 0) then
           writeln(n,x:14,ynew[1]:14,ynew[2]:14,hu:14,hv:14);
    end;
    readln;
end.