Klassiches Runge-Kutta-Verfahren

program RungeKutta(input, output);

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

    var
        x0, y0: double; {Anfangswerte}
        xf:     double; {Endwert}
        h:      double; {Schrittweite}
        x, y, x1, y1, x2, y2, x3, y3: double;
        k1, k2, k3, k4: 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 klassischen Runge-Kutta-Verfahren am Beispiel');
        writeln('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}
        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);
           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
        end; {while}
    end. {Runge Kutta}