(*$B-,U+ Aufgabe 7.29 *)
PROGRAM fehlberg;
CONST nn=10;
TYPE vektor = ARRAY[1..nn] of real;
VAR x,x0,h,xend,xmax, xmin, ymin, ymax,
    man,max,norm,fehler,eps : real ;
    m,n,i,k, fall, xg, yg, xt, yt, oldx, oldy, farbe, keine,miss,gel: integer;
    zeich,misslungen : boolean;
    y,y0,ys,y1,y2,z,k1,k2,k3,k4,k5,k6: vektor;
    tf,aus:text; line : string[120];

 (*$I plotproz *)

 PROCEDURE f(x:real; y:vektor; VAR ys:vektor);
 BEGIN
    CASE fall OF
       (*$I fdgl.pas*)
    END;
 END;

 PROCEDURE schritt;
 BEGIN
    REPEAT
       f(x,y,k1);
       FOR i:=1 TO n DO z[i]:=y[i]+h/4*k1[i];
       f(x+h/4,z,k2);
       FOR i:=1 TO n DO z[i]:=y[i]+h/32*(3*k1[i]+9*k2[i]);
       f(x+3*h/8,z,k3);
       FOR i:=1 TO n DO
         z[i]:=y[i]+h/2197*(1932*k1[i]-7200*k2[i]+7296*k3[i]);
       f(x+12*h/13,z,k4);
       FOR i:=1 TO n DO
         z[i]:=y[i]+h*(439*k1[i]/216-8*k2[i]
                    +3680*k3[i]/513-845*k4[i]/4104);
       f(x+h,z,k5);
       FOR i:=1 TO n DO
       z[i]:=y[i]+h*(-8*k1[i]/27+2*k2[i]-3544*k3[i]/2565
                    +1859*k4[i]/4104-11*k5[i]/40);
       f(x+h/2,z,k6);
       norm:=0;fehler:=0;
       FOR i:=1 TO n DO
       BEGIN
         y1[i]:=y[i]+h*(25*k1[i]/216+1408*k3[i]/2565
                     + 2197*k4[i]/4104-k5[i]/5);
         y2[i]:=y[i]+h*(16*k1[i]/135+6656*k3[i]/12825.0
                   +28561.0*k4[i]/56430.0 -9*k5[i]/50+2*k6[i]/55);
         max:=abs(y1[i]-y2[i]);
         IF max>fehler THEN fehler:=max;
         max:=abs(y2[i]);
         IF max>norm THEN norm:=max;
       END;
       misslungen := fehler > eps*norm;
       IF misslungen THEN
       BEGIN  miss:=miss+1; h := 0.8*h END;
    UNTIL NOT misslungen;
    gel:=gel+1;
 END ;

 BEGIN
    assign(tf,'fdgl.pas'); reset(tf);
    while not eof(tf) do
    begin readln(tf,line); writeln(line)  end;
    writeln('Fall eingeben:'); read(fall);
    writeln('Anzahl Gleichungen ?'); read(n);
    write('Anfangsbedingungen eingeben x='); read(x0);
    FOR i := 1 TO n DO BEGIN write('y[',i,']='); read(y0[i]) END;
    writeln('bis wohin integrieren? xend='); read(xend);
    writeln('EPS =?'); read(eps);
    writeln('Wertetabelle (1) oder Zeichnung (2) ?');
    readln(i); zeich := i=2;
    IF zeich THEN
    BEGIN
      writeln('welche Funktion plotten ?'); read(m);
      writeln('ymin, ymax eingeben'); readln( ymin, ymax);
      xmin:=x0; xmax :=xend;
      plotbegin( xg, yg, xt, yt, farbe, keine);
      achsen(xmin, xmax, ymin, ymax, xg, yg, xt, yt, farbe, keine);
    END
    ELSE
    BEGIN
       writeln('Filenamen fuer die Wertetabelle ?');
       readln(line); assign(aus,line); rewrite(aus);
    END;
    miss:=0; gel:=0; h:=(xend-x0)/100; x:=x0; y:=y0;
    REPEAT
       schritt; y:=y2; x := x+h;
       IF zeich THEN pl(x,y[m],farbe)
       ELSE
       BEGIN
          write(aus,x);
          FOR i:=1 TO n DO write(aus,y[i]);
          writeln(aus)
       END;
       IF fehler=0 THEN h:= 1.2*h
       ELSE h:=0.8*h*exp(ln(eps*norm/fehler)/5) ;
    UNTIL x>=xend;
    IF zeich THEN gotoxy(1,23);
    writeln('misslungen:',miss,' gelungen:',gel,' Schritte');
    IF NOT zeich THEN close(aus)
  END.