数値計算プログラム

input F$;
input h,e;
input n,l;
dim x(n),y(n),f(n);
input x(1),y(1);

for i:=1 to 5;
x:=x(i);
y:=y(i);
c1:=h*val(F$);
x:=x(i)+h/3;
y:=y(i)+c1/3;
c2:=h*val(f$);
x:=x(i)+h*0.4;
y:=y(i)+(6*c2+4*c1)/25;
c3:=h*val(f$);
x:=x(i)+h;
y:=y(i)+(15*c3-12*c2+c1)/4;
c4:=h*val(f$);
x:=x(i)+2*h/3;
y:=y(i)+(8*c4-50*c3+90*c2+6*c1)/81;
c5:=h*val(f$);
x:=x(i)+0.8*h;
y:=y(i)+(8*c4+10*c3+36*c2+6*c1)/75;
c6:=h*val(f$);
x(i+1):=x(i)+h;
y(i+1):=y(i)+(23*c1+125*c3-81*c5+125*c6)/192;
x:=x(i+1);
y:=y(i+1);
f(i+1):=val(f$);
next i;

for i:=6 to n-1;
p1:=1901*f(i)-2774*f(i-1)+2616*f(i-2)-1274*f(i-3)+251*f(i-4);
p:=y(i)+(h/720)*p1;
x(i+1):=x(i)+h;
for j:=1 to l;
x:=x(i+1);
y:=p;
f(i+1):=val(f$);
q1:=251*f(i+1)+646*f(i)-264*f(i-1)+106*f(i-2)-19*f(i-3);
q:=y(i)+h/720*q1;
if abs(p-q) p:=q;
next j;
n:=i;
goto l2;
label l1;
y(i+1):=q;
x:=x(i+1);
y:=y(i+1);
f(i+1):=val(f$);
next i;
label l2;
input i1,i2;
if (i-1)*h for i:=1 to i1*i2/h step i2/h;
print x(i),y(i);
next i;
end;