2. Other methods introduced by the operational calculus.

2.1 The methods using the polynomials of higher degree.
. . The author's operational calculus is able to introduce many methods using the polynomials of higher degree. If four points are equidistant, the predictor of y3 is appended to the procedure[3] of Eqs.(1.2), which are not repeated. Then the procedure[4] corrects y1, y2 and y3 as in Eqs.(2.1) and is repeated twice more. The corrector of y3 is called Simpson's 3/8 rule.
. . If five points are equidistant, the predictor of y4 is appended to the procedure[4] without repeat and the procedure[5] corrects y1, y2, y3 and y4 as in Eqs.(2.1), which are repeated twice more. The predictor of y4 is called Milne's predictor and the last corrector in the procedure[5] is called Newton-Cotes's 5 points rule.

[1]f0=f(x0, y0),y1=y0+hf0
[2]f1=f(x1, y1),y1=y0+h(f0+f1)/2,y2=y0+2hf1
[3]f1=f(x1, y1),f2=f(x2, y2), y1=y0+h(5f0+8f1−f2)/12,y2=y0+2h(f0+4f1+f2)/6,
appendy3=y0+3h(f0+3f2)/4
[4]f1=f(x1, y1),f2=f(x2, y2),f3=f(x3, y3), y1=y0+h(9f0+19f1−5f2+f3)/24
y2=y0+2h(f0+4f1+f2)/6, y3=y0+3h(f0+3f1+3f2+f3)/8
appendy4=y0+4h(2f1−f2+2f3)/3(2.1)
[5]f1=f(x1, y1),f2=f(x2, y2),f3=f(x3, y3),. .f4=f(x4, y4),
y1=y0+h(251f0+646f1−264f2+106f3−19f4)/720, . .y2=y0+2h(29f0+124f1+24f2+4f3−f4)/180
y3=y0+3h(9f0+34f1+24f2+14f3−f4)/80, . .y4=y0+4h(7f0+32f1+12f2+32f3+7f4)/90
. . It is not proper that Milne's corrector is Simpson's 1/3 rule. It should be Newton-Cotes's 5 points rule. It is also not proper that usual P-C methods use the corrector of lower degree than the predictor. Furthermore, usual integral rules of higher order are not complete methods of numerical solution of differential equations. They are parts of the author's methods. Usual P-C methods must be given some starting values except initial values because they are parts of numerical solution.
. . Program 2.1 represents the method using 5 points in Eqs.(2.1) and gives more accurate results than the method using 3 points in usual case.

Program 2.1
10 DIM Y0(6),Y1(6),Y2(6),Y3(6),Y4(6),Y(6):DEFINT E,J
12 DEF FNY1=(5*Y0(J3+1)+8*Y1(J3+1)-Y2(J3+1))*H3+Y0(J3):
. . DEF FNY2=(Y0(J3+1)+4*Y1(J3+1)+Y2(J3+1))*H4+Y0(J3):
. . DEF FNY3=(Y0(J3+1)+3*Y2(J3+1))*H5+Y0(J3)
14 DEF FNY4=(9*Y0(J3+1)+19*Y1(J3+1)-5*Y2(J3+1)+Y3(J3+1))*H6+Y0(J3):
. . DEF FNY6=(Y0(J3+1)+3*Y1(J3+1)+3*Y2(J3+1)+Y3(J3+1))*H8+Y0(J3):
. . DEF FNY7=(2*Y1(J3+1)-Y2(J3+1)+2*Y3(J3+1))*H9+Y0(J3)
16 DEF FNY8=(251*Y0(J3+1)+646*Y1(J3+1)-264*Y2(J3+1)+106*Y3(J3+1)-19*Y4(J3+1))*H10+Y0(J3):
. . DEF FNY9=(29*Y0(J3+1)+124*Y1(J3+1)+24*Y2(J3+1)+4*Y3(J3+1)-Y4(J3+1))*H11+Y0(J3):
. . DEF FNY10=(9*Y0(J3+1)+34*Y1(J3+1)+24*Y2(J3+1)+14*Y3(J3+1)-Y4(J3+1))*H12+Y0(J3)
18 DEF FNY11=(7*Y0(J3+1)+32*Y1(J3+1)+12*Y2(J3+1)+32*Y3(J3+1)+7*Y4(J3+1))*H13+Y0(J3)
20 GOSUB 54:H0=(XE-XB)/N:PRINT"X= ";XB;"-";XE;" N=";N;" H=";H0:N1=ND+1
22 X4=XB:FOR IB=1 TO N:X0=X4:X4=XB+H0*IB:FOR J3=1 TO ND:Y0(J3)=Y(J3):NEXT
24 LOCATE 0,CSRLIN:PRINT"IB=";IB;" X=";X4;:GOSUB 40:X=X0:GOSUB 52:Y0(N1)=F
26 FOR J3=1 TO ND:Y1(J3)=Y0(J3+1)*H+Y0(J3):Y(J3)=Y1(J3):NEXT:GOSUB 42
28 FOR J3=1 TO ND:Y1(J3)=(Y0(J3+1)+Y1(J3+1))*H*.5+Y0(J3):
. . Y2(J3)=Y1(J3+1)*H2+Y0(J3):Y(J3)=Y1(J3):NEXT:GOSUB 42:GOSUB 44
30 FOR J3=1 TO ND:Y1(J3)=FNY1:Y2(J3)=FNY2:Y3(J3)=FNY3:Y(J3)=Y1(J3):NEXT:
. . GOSUB 42:GOSUB 44:GOSUB 46
32 FOR J3=1 TO ND:Y1(J3)=FNY4:Y2(J3)=FNY2:Y3(J3)=FNY6:Y4(J3)=FNY7:
. . Y(J3)=Y1(J3):NEXT:GOSUB 42:GOSUB 44:GOSUB 46:GOSUB 48
34 FOR J3=1 TO ND:Y1(J3)=FNY8:Y2(J3)=FNY9:Y3(J3)=FNY10:Y4(J3)=FNY11:
. . Y(J3)=Y1(J3):NEXT:GOSUB 42:GOSUB 44:GOSUB 46:GOSUB 48
36 FOR J3=ND TO 1 STEP-1:Y1(J3)=FNY8:Y2(J3)=FNY9:Y3(J3)=FNY10:
. . Y4(J3)=FNY11:Y(J3)=Y1(J3):NEXT:GOSUB 42:GOSUB 44:GOSUB 46:
. . GOSUB 48:J3=ND:Y4(J3)=FNY11:Y(ND)=Y4(ND)
38 GOSUB 50:NEXT IB:END
40 HM=X4-X0:H=H0*.25:H2=H0*.5:H1=H0*.75:H3=H/12:H4=H2/6:H5=H1/4:
. . H6=H/24:H8=H1/8:H9=HM/3:H10=H/720:H11=H/90:H12=H1/80:H13=HM/90:
. . X1=X0+H:X2=X0+H2:X3=X0+H1:RETURN
42 X=X1:GOSUB 52:Y1(N1)=F:RETURN
44 X=X2:FOR J3=1 TO ND:Y(J3)=Y2(J3):NEXT:GOSUB 52:Y2(N1)=F:RETURN
46 X=X3:FOR J3=1 TO ND:Y(J3)=Y3(J3):NEXT:GOSUB 52:Y3(N1)=F:RETURN
48 X=X4:FOR J3=1 TO ND:Y(J3)=Y4(J3):NEXT:GOSUB 52:Y4(N1)=F:RETURN
50 X#=X4:GOSUB 56:YY#=Y(1):LOCATE 40,CSRLIN:PRINT USING
. . "Y=##.#######^^^^ True=##.########^^^^";YY#;Y#:RETURN

52 F=100*(SIN(X)-Y(1)):RETURN
54 XB=0:XE=10:N=500:ND=1:Y(1)=0:RETURN:'h=.02
56 Y#=(SIN(X#)-.01#*(COS(X#)-EXP(-100*X#)))/1.0001#:RETURN

2.2 Accuracy of the results by Program 2.1.
. . Table 2.1 shows the results of a stiff differential equation by Program 2.1 and the author's 3 points method. The 5 points method with step interval 0.02 has the same condition h=0.005 for existence of solution as the 3 points method with step interval 0.01 and the results are slightly worse than the results of two steps of the 3 points method. However it is not correct to draw the conclusion that the 5 points method is not better in accuracy because the smaller the difference of solution by one step is, the better its accuracy is, and the difference y(0.02)−y(0) of one step by the 5 points method is 3 times difference y(0.01)−y(0) of the first step by the 3 points method.
Table 2.1 Results of y' =100(sin t −y), y(0)=0
t5 points method3 points
method		Accurate
solution		10n
step=0.01step=0.02step=0.01
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.19
0.20		 3.6785675
1.1352742
2.0495741
3.0177429
4.0055022
5.0001778
5.9970450
6.9943026
7.9912253
8.9874841
9.9828921
1.0977320
1.1970656
1.2962799
1.3953644
1.4943098
1.5931056
1.6917419
1.7902094
1.8884978		 :
1.1305087
:
3.0164769
:
4.9999267
:
6.9942616
:
8.9874819
:
1.0977326
:
1.2962806
:
1.4943106
:
1.6917431
:
1.8884990		 3.6805207
1.1354182
2.0496540
3.0177828
4.0055208
5.0001867
5.9970498
6.9943063
7.9912283
8.9874871
9.9828951
1.0977323
1.1970660
1.2962803
1.3953649
1.4943101
1.5931059
1.6917424
1.7902099
1.8884982		 3.67875972
1.13528838
2.04958203
3.01774709
4.00550414
5.00017871
5.99704594
6.99430354
7.99122609
8.98748430
9.98289234
1.09773202
1.19706571
1.29628001
1.39536462
1.49430988
1.59310573
1.69174212
1.79020948
1.88849783		 -03
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-01
-01
-01
-01
-01
-01
-01
-01
-01
. . If the both methods have the same difference of solution by one step, the 5 points method has better results about two figures than the 3 points method. However the interval h of 5 points method is a half of the 3 points method and this comparison is unfair. It is because the smaller the interval h is, the better the accuracy of solution is, even if the method is same.
. . Therefore it is not able to compare the methods using a different number of points by all the same conditions. However it will be correct to draw the conclusion that the 5 points method has better accuracy than the 3 points method from Table 2.1.
. . If the difference of solution by one step is small, the 5 points method is better than the 3 points method as shown over the range t<5 of the 5 points method with step 0.2 and the 3 points method with step 0.1 in Table 2.2. However the one is slightly worse than the other around t=10. In case of h=0.1, where the condition for existence of solution is satisfied over the range t>6.2, the 5 points method with step 0.4 shows very worse results as if they are Milne's results, whereas the 3 points method with step 0.2 does not show the phenomena as yet.
. . However it is not correct to draw the conclusion that the 3 points method is stable and the 5 points method is unstable, because the 3 points method can not obtain more accurate results even if it starts up at around t=7 with accurate initial value and in this case, there is more strict condition for existence of solution depending on t as mentioned after in the theorem 2.1, which shows the condition y(t0+h)>y(t0)/2 must be satisfied.
Table 2.2 Comparison of results for y' =−ty, y(0)=10
t5 points method3 points methodAccurate
solution		10n
step=0.2step=0.4step=0.1step=0.2
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.4
2.6
2.8
3.0
3.2
3.4
3.6
3.8
4.0
4.2
4.4
4.6
4.8
5.0
5.2
5.4
5.6
5.8
6.0
6.2
6.4
6.6
6.8
7.0
7.2
7.4
7.6
7.8
8.0
8.2
8.4
8.6
8.8
9.0
9.2
9.4
9.6
9.8
10.0		 9.8019867
9.2311630
8.3527012
7.2614899
6.0653062
4.8675218
3.7531109
2.7803729
1.9789866
1.3533529
8.8921607
5.6134742
3.4047434
1.9841090
1.1108987
5.9760120
3.0887051
1.5338017
7.3179612
3.3545773
1.4774466
6.2519166
2.5417918
9.9286321
3.7261707
1.3435525
4.6544174
1.5491348
4.9535777
1.5217796
4.4913250
1.2734457
3.4685890
9.0755092
2.2809560
5.5062621
1.2766163
2.8423704
6.0765867
1.2472022
2.4571213
4.6453811
8.4258602
1.4657175
2.4442777
3.9057119
5.9762683
8.7508212
1.2251042
1.6381593		 :
9.2311630
:
7.2614903
:
4.8675227
:
2.7803733
:
1.3533508
:
5.6133890
:
1.9839504
:
5.9740677
:
1.5320852
:
3.3431742
:
6.1933859
:
9.6927281
:
1.2681434
:
1.3573090
:
1.1355667
(exp)
6.69(-09)
:
2.06(-10)
:
-7.28(-13)
:
3.47(-14)
:
-3.68(-15)
:
6.79(-16)
:
-1.96(-16)
:
8.48(-17)
:
-5.25(-17)
:
4.56(-17)		 9.8019867
9.2311630
8.3527012
7.2614894
6.0653052
4.8675203
3.7531087
2.7803702
1.9789841
1.3533506
8.8921440
5.6134635
3.4047377
1.9841076
1.1109003
5.9760425
3.0887391
1.5338328
7.3182075
3.3547555
1.4775656
6.2526530
2.5422205
9.9309967
3.7274011
1.3441613
4.6572768
1.5504139
4.9590409
1.5240036
4.4999805
1.2766615
3.4800143
9.1143482
2.2935817
5.5455442
1.2883097
2.8756873
6.1674662
1.2709271
2.5164109
4.7872358
8.7506370
1.5368828
2.5935047
4.2051071
6.5509282
9.8055567
1.4101849
1.9485392		 9.8019867
9.2311630
8.3526993
7.2614827
6.0652909
4.8674970
3.7530768
2.7803323
1.9789442
1.3533131
8.8918263
5.6132221
3.4045735
1.9840091
1.1108503
5.9758529
3.0887170
1.5338806
7.3188320
3.3552786
1.4779235
6.2547845
2.5433497
9.9363431
3.7296551
1.3449844
4.6597006
1.5508344
4.9579864
1.5224336
4.4895756
1.2712722
3.4557770
9.0159480
2.2568303
5.4178759
1.2467912
2.7487145
5.8011950
1.1710647
2.2586831
4.1569938
7.2895571
1.2156665
1.9237984
2.8812736
4.0709053
5.4053022
6.71(-21)
7.75(-22)		 9.80198673
9.23116344
8.35270199
7.26149030
6.06530660
4.86752228
3.75311111
2.78037290
1.97898674
1.35335283
8.89216081
5.61347500
3.40474421
1.98410974
1.11089965
5.97602198
3.08871441
1.53380989
7.31802551
3.35462628
1.47748183
6.25214775
2.54193577
9.92949522
3.72665317
1.34381028
4.65571332
1.54975396
4.95639984
1.52299797
4.49634149
1.27540685
3.47589347
9.10145896
2.28973485
5.53459867
1.28533632
2.86797709
6.14838727
1.26641655
2.50622581
4.76528183
8.70539806
1.52797740
2.57675711
4.17501738
6.49931301
9.72094942
1.39694133
1.92874985		 00
00
00
00
00
00
00
00
00
00
-01
-01
-01
-01
-01
-02
-02
-02
-03
-03
-03
-04
-04
-05
-05
-05
-06
-06
-07
-07
-08
-08
-09
-10
-10
-11
-11
-12
-13
-13
-14
-15
-16
-16
-17
-18
-19
-20
-20
-21
. . Table 2.3 shows accurate solutions by the interval h=0.1 and that the above condition is not satisfied where t0 is grater than 6.9. The theorem 2.1 does not show that it is sufficient for the interval h to satisfy the condition. It is desirable ideally that the condition is satisfied over the range of one step. When the interval is h=0.1, the interval of one step is 4h=0.4 in case of 5 points method and the solution becomes 1/10 of the starting value around t=6.0. It causes that the solution has only one significant figure in accuracy and go to the state being said unstable. In case of 3 points method, because the interval of one step is 2h=0.2 and the solution is greater than 1/4 of the starting value, the solution has three significant figures in accuracy. The place where the solution by one step of 5 points method is about 1/4 of the starting value is around t=3.2 and the solution has four significant figures in accuracy. It is said usually that higher order methods are unstable but it is only the cause that the selection of interval is not proper in case of author's method and higher order methods are not unstable. If the step interval is 0.3, that is h=0.15, the 3 points method also shows abnormal results around t=8. Runge-K's method also has similar trouble and all the methods using constant interval has this trouble depending on the solution of differential equation.
Table 2.3
taccurate yexp
6.0
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
6.9
7.0
7.1
7.2		 1.52299797
8.31670729
4.49634149
2.40671954
1.27540685
6.69158609
3.47589347
1.78755544
9.10145896
4.58795947
2.28973485
1.13137839
5.53459867		 -07
-08
-08
-08
-08
-09
-09
-09
-10
-10
-10
-10
-11

. . [Theorem 2.1]. . If the analytical solution of y' =f(y,t) is known, the condition for existence of numerical solution is
Eq.(2.2)
M is the infimum satisfying Eq.(1.7) in proof of Theorem 1.1, if f(y,t) is any unknown function. However if f(y,t) is known, M must be the minimum satisfying Eq.(2.3) because f(y0,t0) is fixed and f(y,t) is iterated.

|f(y(a),a)(t−t0)|≤M|f(y,t)|(t0<a<t)(2.3)
Then M must satisfy |y(t)−y(t0)|=M|y' (t)|.[Q.E.D.]

. . In case of y' =−ty and t>t0, because of y(t)<y(t0), Eq.(2.2) comes to

−y(t)+y(t0)<y(t)∴. .y(t)>y(t0)/2
This condition is not satisfied over the range t0≥6.9 in case of h=0.1 as shown in Table 2.3.
. . In case of y' =100(sin t −y) and t>t0≈0, because of y(t0)<y(t)<sin t, Eq.(2.2) comes to
y(t)−y(t0)<sin t −y(t)∴. .y(t)<{y(t0)+sin t}/2
Although this condition is most severe at y(t0)=0, it is satisfied from the accurate solution in Table 2.1, even if h is 0.01.
. . Theorem 1.4 requires that f(y(a),a) is nearly equal to f(y,t0+h). On the other hand, Theorem 2.1 does not require that condition. Hence it gives the more exact condition for existence of solution. However it requires the solution to be known. Therefore it can not be used to estimate the condition for existence of solution but it is of use in order to research the reason of abnormal results.

2.3 The variable pitch method using 3 points.
. . The stiff differential equation in Eq.(1.16) requires the condition 100h<1 over all the range. However it is not sufficient in order to obtain accurate results in the range where the derivative dt/dy is not negligible. In the neighbourhood of t=0, the condition should be expressed in Eq.(1.17) and it is very severe to obtain accurate results. The stiff differential equation in Eq.(1.18) also requires the severe condition of h<0.000998 and it is more severe in the neighbourhood of t=0 when it is required for the results to be as accurate as the results in other range. However the condition is not obtained because of no explicit relation between t and y' .
. . The author's method is able to obtain the accurate solutions by use of variable pitch technique without calculating the condition of interval. The method is very simple because the procedure[3] in Eqs.(1.2) repeated three times gives the data to be able to detect the proper condition.
. . When these are expressed in r1, r2 and r3, the difference of r1 and r2 is an approximate remainder of Newton's interpolation formula and the difference of r1 and r3 is the remainder corrected once by the procedure[3]. It is meaningless to iterate the procedure[3] more times in the theory of the operational calculus except the case that the solution is lower order polynomial than 4th order or has negligible remainder depending on the interval h. Hence the difference of r2 and r3 may be the error of the remainder approximated first. If it is not negligible, the interval h is not proper. Therefore a proper interval is estimated by the procedure as follows.

(1). .If |r2−r3| is not negligible, the step interval 2h is divided by two.
(2). .If |r2−r3| is negligible, the solution is convergent and then,
(3). .If |r1−r3| is negligible and the number of the step is even, the next two steps are
. . .connected so the step interval comes to twice the last one.

Program 2.2
10 DIM Y0(6),Y1(6),Y2(6),Y(6):DEFINT E,J
12 DEF FNY1=(5*Y0(J3+1)+8*Y1(J3+1)-Y2(J3+1))*H3+
. . Y0(J3):DEF FNY2=(Y0(J3+1)+4*Y1(J3+1)+Y2(J3+1))
. . *H4+Y0(J3)
14 RP=16777216#:RM=2/RP
16 GOSUB 58:H0=(XE-XB)/N:PRINT"X= ";XB;"-";XE;
. . " N=";N;" H=";H0:N1=ND+1
18 EE=1:XL=XB:X=XB:GOSUB 48:FOR IB=1 TO N:ES=1:
. . XI=XL:XL=XB+H0*IB:X0=XI
20 WHILE 2:H1=H0/EE:EB=ES:WHILE EB=<EE:
. . X2=XI+H1*EB:IF EB=EE THEN X2=XL
22 LOCATE 0,CSRLIN:PRINT"IB=";IB;" EB=";EB;
. . " EE=";EE;" X=";X2;:GOSUB 46
24 FOR J3=1 TO ND:Y1(J3)=Y0(J3+1)*H+Y0(J3):
. . Y(J3)=Y1(J3):NEXT:GOSUB 50
26 FOR J3=1 TO ND:Y1(J3)=(Y0(J3+1)+Y1(J3+1))*H*
. . .5+Y0(J3):Y2(J3)=Y1(J3+1)*H2+Y0(J3):
. . Y(J3)=Y1(J3):NEXT:GOSUB 50:GOSUB 52
28 IF ABS(Y2(ND))<1 THEN RA=RP ELSE RA=1
30 FOR J3=1 TO ND:Y1(J3)=FNY1:Y2(J3)=FNY2:
. . Y(J3)=Y1(J3):NEXT:GOSUB 50:GOSUB 52
32 R1=Y(ND)*RA
34 FOR J3=ND TO 1 STEP-1:Y1(J3)=FNY1:Y2(J3)=FNY2
. . :Y(J3)=Y1(J3):NEXT:GOSUB 50:GOSUB 52
36 R2=Y(ND)*RA:J3=ND:Y2(J3)=FNY2:Y(J3)=Y2(J3):
. . R3=Y(ND)*RA
38 RA=RM*ABS(R3):IF ABS(R2-R3)<RA THEN
. . GOSUB 54:X=X2:GOSUB 48:X0=X2:GOSUB 42 ELSE
. . IF EE<&H4000 THEN EE=EE*2:ES=EB*2-1:EB=EE+1
. . ELSE PRINT"EE>&h4000":END
40 EB=EB+1:WEND:IF EB>EE+1 THEN WEND ELSE
. . NEXT IB:END
42 IF ABS(R1-R3)<RA THEN IF EE>1 AND(EB MOD 2=0)
. . THEN EE=EE\2:ES=EB\2+1:IF ES=<EE THEN EB=EE+1
44 RETURN
46 H2=X2-X0:H=H1*.5:H3=H/12:H4=H2/6:X1=X0+H:
. . RETURN
48 GOSUB 56:Y0(N1)=F:FOR J3=1 TO ND:Y0(J3)=Y(J3)
. . :NEXT:RETURN
50 X=X1:GOSUB 56:Y1(N1)=F:RETURN
52 X=X2:FOR J3=1 TO ND:Y(J3)=Y2(J3):NEXT:
. . GOSUB 56:Y2(N1)=F:RETURN
54 X#=X2:GOSUB 60:YY#=Y(1):LOCATE 40,CSRLIN:PRINT
. . USING"Y=##.#######^^^^ True=##.########^^^^";
. . YY#;Y#:RETURN

56 F=100*(SIN(X)-Y(1)):RETURN
58 XB=0:XE=50:N=500:ND=1:Y(1)=0:RETURN
60 Y#=(SIN(X#)-.01#*(COS(X#)-EXP(-100*X#)))/
. . 1.0001#:RETURN
 [Note] Y(k+1)=y(k), ND=n-th order.
N=the number of equidistances between XB to XE.

. . In above procedure (3), there are two reasons why the next two steps are connected when the number of the step is even. One of them is exactly to recover the interval before. The other is to avoid bistable repeat of dividing and connecting an interval.
. . In case of n-th order differential equation, r1, r2 and r3 are the values y(n−1) iterated three times. Other values from y(n−2) to y are the parameters of the right side function of the first order differential equation as to y(n−1). So these are corrected twice times.
. . Program 2.2 represents the variable pitch method in BASIC. Line 38 is the above procedure (1) and (2), and line 42 is the procedure (3). Line 28 multiplies r1, r2 and r3 by 224 in case of |y2|<1 in the procedure[2] of Eqs.(1.2) for the purpose of avoiding underflow of exponent by these differences. "WHILE 2" in line 20 is the infinite repeat which ends with "ELSE" in line 40.
. . Line 56, 58 are the subroutines for solving the stiff differential equation in Eq.(1.16) and Table 2.4 represents a part of the results. The basic interval of a step is equidistant and it does not come to grater one than 2h=0.1 because it is required to divide the total range [0,50] into 500. If the results are required only at these equidistant points, line 54 may be appended "IF EB<EE THEN RETURN ELSE " to the head.

2.4 Results by the variable pitch method.
Table 2.4 Solution of y' =100(sin t −y), y(0)=0
by the variable pitch method.
IBEBEEt (10n)yAccurate y10n
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2		 1
2
3
14
8
24
13
32
17
15
16
1
2
5
10
6
16		 512
512
512
512
256
256
128
128
64
16
16
16
16
32
32
16
16		 1.95312(-4)
3.90625(-4)
5.85937(-4)
2.73438(-3)
0.003125
0.009375
0.0101563
0.025
0.0265625
0.09375
0.1
0.10625
0.1125
0.115625
0.13125
0.1375
0.2		 1.8949914
7.5310159
1.6835715
3.4197359
4.4115592
3.2910295
3.7779736
1.5819790
1.7263290
8.3649129
8.9874841
9.6097231
1.0231596
1.0542386
1.2094744
1.2714882
1.8884976		 1.89499145
7.53101618
1.68357167
3.41973592
4.41155928
3.29102951
3.77797369
1.58197891
1.72632887
8.36491276
8.98748430
9.60972308
1.02315963
1.05423863
1.20947437
1.27148822
1.88849783		 -6
-6
-5
-4
-4
-3
-3
-2
-2
-2
-2
-2
-1
-1
-1
-1
-1
. . Although line 58 in Program 2.2 requires to solve the differential equation with the step interval 2h=0.1, the solution does not converge with the step interval because it does not satisfy the condition in Eq.(1.17).The solution converges first at t=0.1/512 when the step interval has been divided into 512. The step interval is connected and comes to 0.1/256 after 14th step, which is 7th step of division into 256. The step interval is connected and comes to 0.1/128 after 24th step, which is 12th step of division into 128. Thus the step interval comes to division into 16 around t=0.1 which is the first point of the basic division. The step interval goes on with division into 16 or 32 since then, because it sufficiently satisfies the condition h<0.01 but the division into eight is not sufficient for obtaining the results with single precision.
. . The differential equation in Eq.(2.4) is also stiff and requires the condition h<0.000998. However it is not sufficient for obtaining the accurate results with single precision in the neighbourhood of t=0 and the sufficient condition is not obtained. Program 2.2 can obtain the accurate results without use of the condition by use of following subroutine, which requires the results by the interval 0.1. Table 2.5 and 2.6 show the results.

y" +1001y' +1000y=0,y1(0)=1, y1' (0)=998,y2(0)=0, y2' (0)=−999(2.4)
56 F=-1001*Y(2)-1000*Y(1):RETURN
58 XB=0:XE=5:N=50:ND=2:Y(2)=998:Y(1)=1:RETURN
60 Y#=2*EXP(-X#)-EXP(-1000*X#):RETURN		:' For y2
:' Y(2)=-999:Y(1)=0
:' Y#=-EXP(-X#)+EXP(-1000*X#)

. . Although line 58 in the above subroutine requires to solve the differential equation with the step interval 2h=0.1, the solution does not converge with the step interval. The solution converges first at t=0.1/1024 when the step interval has been divided into 1024. The 21st step does not converge and is divided by 2, where the first step is the 41st step of 2048 steps. The 127th step does not converge and is divided by 2, where the first step is the 253rd step of 4096 steps. The step intervals are connected and come to two times the last interval after 260th step of 4096 steps, 148th step of 2048 steps, 104th step of 1024 steps and so on. Thus the step interval comes to about 0.1/64 around t=0.1 and afterwards. Table 2.6 shows that the solution y2 is obtained in a similar way.
Table 2.5 Results of y1(0)=1, y1' (0)=998
IBEBEEty1Accurate10n
1
1
1
1
1
1
1
1
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19		 1
41
253
131
75
53
35
23
55
64
32
128
16
64
64
64
64
64
64
32
64
32
64
128
64
64
64
64		 1024
2048
4096
2048
1024
512
256
128
64
64
32
128
16
64
64
64
64
64
64
32
64
32
64
128
64
64
64
64		 9.765(-5)
2.001(-3)
6.176(-3)
6.396(-3)
7.324(-3)
0.0103516
0.0136719
0.0179687
0.0859375
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9		 1.0928441
1.8609290
1.9856073
1.9855803
1.9847459
1.9793717
1.9728411
1.9643834
1.8353033
1.8096751
1.6374615
1.4816370
1.3406407
1.2130615
1.0976236
9.9317110
8.9865828
8.1313980
7.3575932
6.6574252
6.0238856
5.4506361
4.9319416
4.4626060
4.0379316
3.6536717
3.3059785
2.9913738		 1.09284408
1.86092888
1.98560741
1.98558045
1.98474570
1.97937172
1.97284117
1.96438344
1.83530317
1.80967483
1.63746150
1.48163642
1.34064008
1.21306132
1.09762325
9.93170619
8.98657918
8.13139290
7.35758882
6.65742152
6.02388395
5.45063547
4.93193940
4.46260320
4.03793026
3.65367031
3.30597753
2.99137246		 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
Table 2.6 Results of y2(0)=0, y2' (0)=−999
IBEBEEty2Accurate10n
1
1
1
1
1
1
1
1
1
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19		 1
43
281
146
81
57
36
24
15
29
32
32
128
16
64
32
64
64
128
64
32
64
32
64
128
64
64
64
64		 1024
2048
4096
2048
1024
512
256
128
64
32
32
32
128
16
64
32
64
64
128
64
32
64
32
64
128
64
64
64
64		 9.765(-5)
2.099(-3)
6.860(-3)
7.128(-3)
7.910(-3)
0.0111328
0.0140625
0.01875
0.0234375
0.090625
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9		 -9.2941724
-8.7539840
-9.9211478
-9.9209505
-9.9175429
-9.8891443
-9.8603523
-9.8142475
-9.7683501
-9.1336042
-9.0483761
-8.1873089
-7.4081856
-6.7032045
-6.0653090
-5.4881197
-4.9658567
-4.4932920
-4.0656993
-3.6787966
-3.3287126
-3.0119428
-2.7253181
-2.4659708
-2.2313030
-2.0189658
-1.8268359
-1.6529892
-1.4956869		 -9.29417319
-8.75398332
-9.92114582
-9.92094848
-9.91754049
-9.88914303
-9.86035133
-9.81424680
-9.76835025
-9.13360154
-9.04837417
-8.18730751
-7.40818212
-6.70320042
-6.06530660
-5.48811623
-4.96585310
-4.49328959
-4.06569645
-3.67879441
-3.32871076
-3.01194198
-2.72531774
-2.46596970
-2.23130160
-2.01896513
-1.82683515
-1.65298876
-1.49568623		 -02
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01

Table 2.7 y' =−ty, y(0)=10, basic H=0.1
IBEBEEtyAccurate10n
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
25
30
35
40
50
60
70
80
90
100
110
120
130		 1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
4
4
4
8
8
8
8
8
16
16
16
16
16		 1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
4
4
4
8
8
8
8
8
16
16
16
16
16		 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0
2.5
3.0
3.5
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0		 9.9501247
9.8019867
9.5599747
9.2311630
8.8249683
8.3527012
7.8270445
7.2614894
6.6697669
6.0653052
5.4607425
4.8675208
4.2955723
3.7531102
3.2465239
2.7803721
2.3574600
1.9789861
1.6447442
1.3533524
4.3936917
1.1108992
2.1874905
3.3546255
3.7266538
1.5229992
2.2897394
1.2664200
2.5767636
1.9287552
5.3111059
5.3802019
2.0050166		 9.95012479
9.80198673
9.55997478
9.23116344
8.82496903
8.35270199
7.82704545
7.26149030
6.66976789
6.06530660
5.46074412
4.86752228
4.29557318
3.75311111
3.24652467
2.78037290
2.35746057
1.97898674
1.64474464
1.35335283
4.39369336
1.11089965
2.18749112
3.35462628
3.72665317
1.52299797
2.28973485
1.26641655
2.57675711
1.92874985
5.31109225
5.38018616
2.00500878		 0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-01
-01
-02
-03
-05
-07
-10
-13
-17
-21
-26
-31
-36
. . Table 2.7 shows the solution of y' =−ty, y(0)=10 with intervals of 2h=0.1. Although the interval h sufficiently satisfies the condition for existence of solution, the result at t=1.2 has the difference 2×10−7 from the accurate solution as shown in the column of 3 points method with step=0.1 in Table 2.2. So the result at t=1.2 is obtained by 2 steps with a half of interval. Around t=10, the results is obtained by 16 steps with the interval 0.1/16 and have the accuracy with 6 significant figures.
. . So this method obtains stably accurate results without estimating the condition of interval h. However it is impossible as usual method to obtain the solutions over the range up to far point from zero, when the solutions decrease more rapidly or increase more slowly than the impulsive response.

2.5 The variable pitch method using 5 points.
. . The author's 5 points method is represented in Eqs.(2.1) and the value y4 of the procedure [5] is corrected three times. Denoting these values as r1, r2, r3, the convergence of the solution and the proper interval h are determined in the same way as the variable pitch method using 3 points.

[1] If |r2−r3| is not negligible, the step interval 4h is divided by two.
[2] If |r2−r3| is negligible, the solution is convergent and then,
[3] If |r1−r3| is negligible and the number of the step is even, the next two
. . steps are connected so the step interval comes to twice the last one.

. . Hence the program of the variable pitch method is constructed by appending the instructions for iteration in Program 2.2 to Program 2.1. The program is represented in Program 2.3. The condition for convergence is that the difference of iterated values is zero to the least two bits and the condition for connecting two step intervals is that it is zero to the least bit. These values have been determined by many trial.
. . Line 66, 68 are the subroutines for solving the given differential equation. The basic interval of a step is equidistant and it does not come to grater one than 4h=0.1 because it is required to divide the total range [0, 13] into 130. If the results are required only at these equidistant points, line 64 may be appended "IF EB<EE THEN RETURN ELSE " to the head.

Program 2.3
10 DIM Y0(6),Y1(6),Y2(6),Y3(6),Y4(6),Y(6):DEFINT E,J
12 DEF FNY1=(5*Y0(J3+1)+8*Y1(J3+1)-Y2(J3+1))*H3+Y0(J3): DEF FNY2=(Y0(J3+1)+4*
. . Y1(J3+1)+Y2(J3+1))*H4+Y0(J3): DEF FNY3=(Y0(J3+1)+3*Y2(J3+1))*H5+Y0(J3)
14 DEF FNY4=(9*Y0(J3+1)+19*Y1(J3+1)-5*Y2(J3+1)+Y3(J3+1))*H6+Y0(J3):
. . DEF FNY6=(Y0(J3+1)+3*Y1(J3+1)+3*Y2(J3+1)+Y3(J3+1))*H8+Y0(J3):
. . DEF FNY7=(2*Y1(J3+1)-Y2(J3+1)+2*Y3(J3+1))*H9+Y0(J3):'fny5=fny2
16 DEF FNY8=(251*Y0(J3+1)+646*Y1(J3+1)-264*Y2(J3+1)+106*Y3(J3+1)-19*Y4(J3+1))*
. . H10+Y0(J3): DEF FNY9=(29*Y0(J3+1)+124*Y1(J3+1)+24*Y2(J3+1)+4*Y3(J3+1)-
. . Y4(J3+1))*H11+Y0(J3): DEF FNY10=(9*Y0(J3+1)+34*Y1(J3+1)+24*Y2(J3+1)+14*
. . Y3(J3+1)-Y4(J3+1))*H12+Y0(J3)
18 DEF FNY11=(7*Y0(J3+1)+32*Y1(J3+1)+12*Y2(J3+1)+32*Y3(J3+1)+7*Y4(J3+1))*H13+
. . Y0(J3):RP=16777216#:RM=4/RP
20 GOSUB 68:H0=(XE-XB)/N:PRINT"X= ";XB;"-";XE;" N=";N;" H=";H0:N1=ND+1
22 EE=1:XL=XB:X=XB:GOSUB 54:FOR IB=1 TO N:ES=1:XI=XL:XL=XB+H0*IB:X0=XI
24 WHILE 2:HE=H0/EE:EB=ES:WHILE EB=<EE:X4=XI+HE*EB:IF EB=EE THEN X4=XL
26 LOCATE 0,CSRLIN:PRINT"IB=";IB;" EB=";EB;" EE=";EE;" X=";X4;:GOSUB 52
28 FOR J3=1 TO ND:Y1(J3)=Y0(J3+1)*H+Y0(J3):Y(J3)=Y1(J3):NEXT:GOSUB 56
30 FOR J3=1 TO ND:Y1(J3)=(Y0(J3+1)+Y1(J3+1))*H*.5+Y0(J3):Y2(J3)=Y1(J3+1)*H2+
. . Y0(J3):Y(J3)=Y1(J3):NEXT:GOSUB 56:GOSUB 58
32 FOR J3=1 TO ND:Y1(J3)=FNY1:Y2(J3)=FNY2:Y3(J3)=FNY3:Y(J3)=Y1(J3):NEXT:
. . GOSUB 56:GOSUB 58:GOSUB 60
34 FOR J3=1 TO ND:Y1(J3)=FNY4:Y2(J3)=FNY2:Y3(J3)=FNY6:Y4(J3)=FNY7:Y(J3)=Y1(J3)
. . :NEXT:GOSUB 56:GOSUB 58:GOSUB 60:GOSUB 62
36 IF ABS(Y4(ND))<1 THEN RA=RP ELSE RA=1
38 FOR J3=1 TO ND:Y1(J3)=FNY8:Y2(J3)=FNY9:Y3(J3)=FNY10:Y4(J3)=FNY11:Y(J3)=
. . Y1(J3):NEXT:GOSUB 56:GOSUB 58:GOSUB 60:GOSUB 62:R1=Y(ND)*RA
40 FOR J3=ND TO 1 STEP-1:Y1(J3)=FNY8:Y2(J3)=FNY9:Y3(J3)=FNY10:Y4(J3)=FNY11:
. . Y(J3)=Y1(J3):NEXT:GOSUB 56:GOSUB 58:GOSUB 60:GOSUB 62:R2=Y(ND)*RA
42 J3=ND:Y4(J3)=FNY11:Y(J3)=Y4(J3):R3=Y(ND)*RA:RA=RM*ABS(R3)
44 IF ABS(R2-R3)<RA THEN GOSUB 64:X=X4:GOSUB 54:X0=X4:GOSUB 48 ELSE
. . IF EE<&H4000 THEN EE=EE*2:ES=EB*2-1:EB=EE+1 ELSE PRINT"EE>&h4000":END
46 EB=EB+1:WEND:IF EB>EE+1 THEN WEND ELSE NEXT IB:END
48 IF ABS(R1-R3)*2<RA THEN IF EE>1 AND(EB MOD 2=0)THEN EE=EE\2:ES=EB\2+1:
. . IF ES=<EE THEN EB=EE+1
50 RETURN
52 HM=X4-X0:H=HE*.25:H2=HE*.5:H1=HE*.75:H3=H/12:H4=H2/6:H5=H1/4:H6=H/24:
. . H8=H1/8:H9=HM/3:H10=H/720:H11=H/90:H12=H1/80:H13=HM/90:X1=X0+H:X2=X0+H2:
. . X3=X0+H1:RETURN
54 GOSUB 66:Y0(N1)=F:FOR J3=1 TO ND:Y0(J3)=Y(J3):NEXT:RETURN
56 X=X1:GOSUB 66:Y1(N1)=F:RETURN
58 X=X2:FOR J3=1 TO ND:Y(J3)=Y2(J3):NEXT:GOSUB 66:Y2(N1)=F:RETURN
60 X=X3:FOR J3=1 TO ND:Y(J3)=Y3(J3):NEXT:GOSUB 66:Y3(N1)=F:RETURN
62 X=X4:FOR J3=1 TO ND:Y(J3)=Y4(J3):NEXT:GOSUB 66:Y4(N1)=F:RETURN
64 X#=X4:GOSUB 70:YY#=Y(1):LOCATE 40,CSRLIN:PRINT USING
. . "Y=##.#######^^^^ True=##.########^^^^";YY#;Y#:RETURN

66 F=-X*Y(1):RETURN
68 XB=0:XE=13:N=130:ND=1:Y(1)=10:RETURN
70 Y#=10*EXP(-X#*X#*.5):RETURN
[Note] Y(k+1)=y(k), ND=n-th order, N=the number of equidistances between XB to XE.


2.6 The results of the variable pitch method using 5 points.
. . The solution of the differential equation given above subroutines is represented in Table 2.8, which has the results of the 3 points method in order to compare the both accuracy. The values of the solutions have no difference but the 5 points method is favourable at less division of the given interval. The 5 points method has the step interval of 4h, where the 3 points method has the step interval of 2h. So if the step interval is two times that of the 3 points method, the both methods have the same value h. However, the table shows that the step interval of the 5 points method comes to four times that of the 3 points method around grater t than t=3.0. It shows that the 5 points method has the better accuracy than the 3 points method.

Table 2.8. . y' = −ty, y(0)=10, basic H=0.1
5 points method3 points method
IBEBEEtyAccurate yexpIBEBEEty
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
25
30
35
40
50
60
70
80
90
100
110
120
130		 1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
4
4
4
4
4
4
8		 1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
4
4
4
4
4
4
8		 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0
2.5
3.0
3.5
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0		 9.9501247
9.8019867
9.5599747
9.2311630
8.8249683
8.3527012
7.8270445
7.2614894
6.6697669
6.0653057
5.4607434
4.8675218
4.2955728
3.7531106
3.2465243
2.7803726
2.3574603
1.9789865
1.6447445
1.3533527
4.3936929
1.1108997
2.1874912
3.3546258
3.7266531
1.5229985
2.2897369
1.2664183
2.5767609
1.9287527
5.3111016
5.3801972
2.0050135		 9.95012479
9.80198673
9.55997478
9.23116344
8.82496903
8.35270199
7.82704545
7.26149030
6.66976789
6.06530660
5.46074412
4.86752228
4.29557318
3.75311111
3.24652467
2.78037290
2.35746057
1.97898674
1.64474464
1.35335283
4.39369336
1.11089965
2.18749112
3.35462628
3.72665317
1.52299797
2.28973485
1.26641655
2.57675711
1.92874985
5.31109225
5.38018616
2.00500878		 +00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
-01
-01
-02
-03
-05
-07
-10
-13
-17
-21
-26
-31
-36		 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
25
30
35
40
50
60
70
80
90
100
110
120
130		 1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
4
4
4
8
8
8
8
8
16
16
16
16
16		 1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
4
4
4
8
8
8
8
8
16
16
16
16
16		 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0
2.5
3.0
3.5
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0		 9.9501247
9.8019867
9.5599747
9.2311630
8.8249683
8.3527012
7.8270445
7.2614894
6.6697669
6.0653052
5.4607425
4.8675208
4.2955723
3.7531102
3.2465239
2.7803721
2.3574600
1.9789861
1.6447442
1.3533524
4.3936917
1.1108992
2.1874905
3.3546255
3.7266538
1.5229992
2.2897394
1.2664200
2.5767636
1.9287552
5.3111059
5.3802019
2.0050166

. . Table 2.9 shows the solutions of the stiff differential equation which is represented in following subroutines. In this case, the 5 points method has eight times the step interval of the 3 points method around t=0. This shows that the stiff differential equations must be solved by the higher order method, if we would like to solve them by using grater interval h, which must satisfy the condition for existence of solutions mentioned in Section 1.5.
. . In case of the fixed pitch method mentioned in Section 2.2, the 5 points method had slightly worse results than the 3 points method. The reason is that the value h is large and near to the critical condition for existence of solutions. Table 2.9 shows that if the value h is smaller and about 0.1/64, the 5 points method has better results than the 3 points method.

66 F=100*(SIN(X)-Y(1)):RETURN
68 XB=0:XE=50:N=500:ND=1:Y(1)=0:RETURN
70 Y#=(SIN(X#)-.01#*(COS(X#)-EXP(-100*X#)))/1.0001#:RETURN

Table 2.9 . .Solution of . .y' =100(sin t −y), y(0)=0, basic H=0.1.
5 points method3 points method
IBEBEEtyAccurate yexpIBEBEEt(10n)yexp
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
3
4
5
6
7
8
9
10		 1
2
3
4
3
4
5
6
7
8
9
10
6
7
8
9
10
6
7
8
1
2
3
4
5
6
7
8
8
8
8
8
8
8
8
8		 64
64
64
64
32
32
32
32
32
32
32
32
16
16
16
16
16
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8		 0.0015625
0.003125
0.0046875
0.00625
0.009375
0.0125
0.015625
0.01875
0.021875
0.025
0.028125
0.03125
0.0375
0.04375
0.05
0.05625
0.0625
0.075
0.0875
0.1
0.1125
0.125
0.1375
0.15
0.1625
0.175
0.1875
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0		 1.1595325
4.4115595
9.4533840
1.6026088
3.2910297
5.3649675
7.7209268
1.0283181
1.2996316
1.5819788
1.8723922
2.1686995
2.7730646
3.3868104
4.0055037
4.6267595
5.2492894
6.4956851
7.7420481
8.9874841
1.0231596
1.1474132
1.2714882
1.3953647
1.5190233
1.6424446
1.7656091
1.8884978
2.8593826
3.8016972
4.7060263
5.5633348
6.3650554
7.1031797
7.7703309
8.3598429		 1.15953252
4.41155928
9.45338329
1.60260869
3.29102951
5.36496733
7.72092668
1.02831811
1.29963169
1.58197891
1.87239224
2.16869970
2.77306478
3.38681073
4.00550414
4.62675988
5.24928980
6.49568572
7.74204870
8.98748430
1.02315963
1.14741320
1.27148822
1.39536477
1.51902335
1.64244460
1.76560908
1.88849783
2.85938259
3.80169721
4.70602653
5.56333504
6.36505606
7.10318000
7.77033129
8.35984363		 -04
-04
-04
-03
-03
-03
-03
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01		 1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
3
4
5
6
7
8
9
10		 1
2
3
14
8
24
13
24
28
32
18
20
24
28
16
18
20
12
14
16
2
8
6
8
10
12
14
16
32
32
32
16
16
32
32
16		 512
512
512
512
256
256
128
128
128
128
64
64
64
64
32
32
32
16
16
16
16
32
16
16
16
16
16
16
32
32
32
16
16
32
32
16		 1.953(-4)
3.906(-4)
5.859(-4)
2.734(-3)
0.003125
0.009375
0.0101563
0.01875
0.021875
0.025
0.028125
0.03125
0.0375
0.04375
0.05
0.05625
0.0625
0.075
0.0875
0.1
0.1125
0.125
0.1375
0.15
0.1625
0.175
0.1875
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0		 1.8949914
7.5310159
1.6835715
3.4197359
4.4115592
3.2910295
3.7779736
1.0283181
1.2996317
1.5819790
1.8723924
2.1686997
2.7730649
3.3868108
4.0055040
4.6267599
5.2492898
6.4956859
7.7420488
8.9874841
1.0231596
1.1474132
1.2714882
1.3953647
1.5190233
1.6424444
1.7656089
1.8884976
2.8593826
3.8016969
4.7060263
5.5633348
6.3650554
7.1031797
7.7703309
8.3598435		 -06
-06
-05
-04
-04
-03
-03
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-02
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01

. . Table 2.10 shows the solutions of the stiff differential equation which is represented in following subroutines. Also in this case, the 5 points method has eight times the step interval of the 3 points method around t=0. This also shows that the stiff differential equations must be solved by the higher order method, if we would like to solve them by using grater interval h.
. . These examples show that it was not correct that the numerical solution would always use a constant order method with a constant interval. It should be the method using variable order polynomials with variable intervals. The author's operational calculus is able to introduce such methods, which will be mentioned afterwards.

66 F=-1001*Y(2)-1000*Y(1):RETURN
68 XB=0:XE=5:N=50:ND=2:Y(1)=1:Y(2)=998:RETURN
70 Y#=2*EXP(-X#)-EXP(-1000*X#):RETURN

Table 2.10 . .y" +1001y' +1000y=0, y(0)=1, y' (0)=998, basic H=0.1.
5 points method3 points method
IBEBEEt(10n)yAccurate yexpIBEBEEt(10n)y
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20		 1
2
3
4
5
36
19
30
16
22
12
16
9
14
29
32
32
32
32
64
64
64
32
64
64
64
64
32
32
32
32
32
32
32
32
64		 512
512
512
512
512
512
256
256
128
128
64
64
32
32
64
64
32
32
32
64
64
64
32
64
64
64
64
32
32
32
32
32
32
32
32
64		 1.953(-4)
3.906(-4)
5.859(-4)
7.812(-4)
9.765(-4)
7.031(-3)
7.421(-3)
0.0117188
0.0125
0.0171875
0.01875
0.025
0.028125
0.04375
0.0453125
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0		 1.1770319
1.3225851
1.4422448
1.5406048
1.6214445
1.9851032
1.9846134
1.9766912
1.9751520
1.9659189
1.9628496
1.9506202
1.9445339
1.9143869
1.9113979
1.9024594
1.8096752
1.6374617
1.4816371
1.3406409
1.2130620
1.0976241
9.9317122
8.9865839
8.1313986
7.3575944
6.6574264
6.0238880
5.4506403
4.9319440
4.4626072
4.0379336
3.6536735
3.3059800
2.9913747
2.7067077		 1.17703185
1.32258506
1.44224466
1.54060475
1.62144438
1.98510300
1.98461317
1.97669115
1.97515187
1.96591869
1.96284937
1.95061982
1.94453365
1.91438645
1.91139756
1.90245885
1.80967483
1.63746150
1.48163642
1.34064008
1.21306132
1.09762325
9.93170619
8.98657918
8.13139290
7.35758882
6.65742152
6.02388395
5.45063547
4.93193940
4.46260320
4.03793026
3.65367031
3.30597753
2.99137246
2.70670566		 +00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
+00
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01
-01		 1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20		 1
41
253
131
148
75
76
60
64
44
24
16
18
14
15
32
64
32
128
16
64
64
64
64
64
64
32
64
32
64
128
64
64
64
64
64		 1024
2048
4096
2048
2048
1024
1024
512
512
256
128
64
64
32
32
64
64
32
128
16
64
64
64
64
64
64
32
64
32
64
128
64
64
64
64
64		 9.765(-5)
2.001(-3)
6.176(-3)
6.396(-3)
7.226(-3)
7.324(-3)
7.421(-3)
0.0117188
0.0125
0.0171875
0.01875
0.025
0.028125
0.04375
0.046875
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0		 1.0928441
1.8609290
1.9856073
1.9855803
1.9848721
1.9847459
1.9846133
1.9766912
1.9751519
1.9659187
1.9628493
1.9506197
1.9445336
1.9143865
1.9084134
1.9024590
1.8096751
1.6374615
1.4816370
1.3406407
1.2130615
1.0976236
9.9317110
8.9865828
8.1313980
7.3575932
6.6574252
6.0238856
5.4506361
4.9319416
4.4626060
4.0379316
3.6536717
3.3059785
2.9913738
2.7067062

TO 2.7 The definite integral by the numerical solution of differential equations.