Verfahren von Heun

program Heun(input, output);

        {Integration der Gleichung y'(x)=-y(x) mittels des Verfahren von
        Heun}

    var
        x0, y0: double; {Anfangswerte}
        h:      double; {Schrittweite}
        eps:    double; {Genauigkeit}
        xf:     double; {Endwert}
        s, s1, x1, y1, del, yneu, x, y: double;
        i, imax:                        integer;

    function fxy(x, y: double): double; begin
        fxy := -y
    end; {fxy}

    procedure init; begin
        writeln('Darstellung der Loesung einer expliziten Differentialgleichung');
        writeln('mittels des Verfahren von Heun am Beispiel y`(x)=-y(x)');
        writeln;
        writeln('Geben Sie bitte die Anfangswerte ein');
        readln(x0, y0);
        writeln('Geben Sie nun an, bis zu welchem x-Wert y(x) berechnet werden soll.');
        readln(xf);
        writeln('Und in welcher Schrittweite sollen die x-Werte liegen?');
        readln(h);
        while (xf - x0)*h < 0 do begin
            writeln('Von ', x0:1:3,' aus koennen Sie ', xf:1:3,' in Schritten von ', h:1:3,' nicht erreichen.');
            writeln('Geben Sie noch einmal den Endwert und die Schrittweite ein.');
            readln(xf, h)
        end; {do}
        writeln('Bis zu welcher Genauigkeit soll der Korrektorschritt durchgefuehrt werden?');
        readln(eps);
        while eps <= 0 do begin
           writeln('Nur positive Genauigkeiten machen Sinn. Versuchen Sie es noch einmal.');
           readln(eps)
        end; {do}
        writeln('Wieviele Iterationen sollen im Korrektorschritt maximal ausgefuehrt werden?');
        readln(imax);
        x := x0;
        y := y0;
        writeln('x       y')
    end; {init}

    procedure result; begin
        writeln( x:1:5,' ', y:1:6)
    end; {result}

    begin {Heun}
        init;
        while x < xf do begin
            s   := fxy(x, y);                      {Praediktor}
            x1  := x + h;
            y1  := y + s*h;
            i   := 0;
            del := 2.0*eps;
            while ((del > eps) and (i < imax)) do begin  {Korrektor}
                 i    := i+1;
                 s1   := fxy(x1, y1);
                 yneu := y + h*(s + s1)/2.0;
                 del  := abs((yneu - y)/yneu);
                 y1   := yneu
                 end; {while}
            if (del > eps) then writeln('der Korrektorschritt konvergiert nicht.');
            x := x1;
            y := yneu;
            result
        end; {while}
    end. {Heun}