Runge-Kutta-Verfahren mit Schrittweitenanpassung

program RungeKuttaSchrittweiteangepasst(input, output);

        {Integration der Gleichung y'(x) = 4*Exp(0.8x) - 0.5y(x) mittels
        des Runge-Kutta Verfahren 4.Ordnung mit Schrittweiteanpassung}

    var
        x0, y0: double; {Anfangswerte}
        xf:     double; {Endwert}
        h:      double; {Schrittweite}
        x, y, x1, y1, x2, y2, x3, y3: double;
        k1, k2, k3, k4, q: double;

    function fxy(x, y: double): double; begin
        fxy := 4.0*Exp(0.8*x) - 0.5*y
    end; {fxy}

    procedure init; begin
        writeln('Darstellung der Loesung einer expliziten Differentialgleichung');
        writeln('mittels des Runge-Kutta-Verfahren mit Schrittweiteanpassung am');
        writeln('Beispiel y`(x) = 4*Exp(0.8x) - 0.5y(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; {while}
        x := x0;
        y := y0;
        writeln('x       y')
    end; {init}

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

    begin {Runge-Kutta-Schr}
        init;
        while x < xf do begin
           k1 := fxy(x, y);
           x1 := x + h/2.0;
           y1 := y + k1*h/2.0;
           k2 := fxy(x1, y1);
           x2 := x1;
           y2 := y + k2*h/2.0;
           k3 := fxy(x2, y2);
           q  := Abs((k3 - k2)/(k2 - k1));
           if q < 0.1 then begin
               x3 := x + h;
               y3 := y + k3*h;
               k4 := fxy(x3, y3);
               y  := y + (1.0/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4)*h;
               x  := x + h;
               result;
               if q < 0.025 then begin
                   h :=2.0*h
               end; {if}
           end else begin {if}
                   h :=h/2.0
           end {if}
        end {while}
    end. {Runge-Kutta-Schr}