'********************************************************************
'*    Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the   *
'*    Adams-Bashforth Method                                        *
'* ---------------------------------------------------------------- *
'* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
'*             Basic et en Pascal By Claude Nowakowski, Edition du  *
'*             P.S.I., 1984" [4].                                   *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'* (Solve Y' = -Y + X/((1+X)*(1+X))  with Y(0)=1 for X0=0 up to     *
'*  X1=1.0, exact solution is Y = 1 / (1+X).                        *
'*                                                                  *
'*      X           Y        Y exact      Error                     *
'* -----------------------------------------------                  *
'*   0.000000    1.000000    1.000000   0                           *
'*   0.000000    1.000000    1.000000   0.000000                    *
'*   0.050000    0.952381    0.952381   0.000000                    *
'*   0.100000    0.909091    0.909091   0.000000                    *
'*   0.150000    0.869525    0.869565   0.000040                    *
'*   0.200000    0.833265    0.833333   0.000068                    *    
'*   0.250000    0.799910    0.800000   0.000090                    *
'*   0.300000    0.769125    0.769231   0.000106                    *
'*   0.350000    0.740623    0.740741   0.000117                    *
'*   0.400000    0.714160    0.714286   0.000125                    *
'*   0.450000    0.689525    0.689655   0.000131                    *
'*   0.500000    0.666533    0.666667   0.000134                    *
'*   0.550000    0.645026    0.645161   0.000135                    *
'*   0.600000    0.624865    0.625000   0.000135                    *
'*   0.650000    0.605926    0.606061   0.000134                    *
'*   0.700000    0.588103    0.588235   0.000133                    *
'*   0.750000    0.571298    0.571429   0.000131                    *
'*   0.800000    0.555428    0.555556   0.000128                    *
'*   0.850000    0.540416    0.540541   0.000125                    *
'*   0.900000    0.526194    0.526316   0.000121                    *
'*   0.950000    0.512703    0.512821   0.000118                    *
'*   1.000000    0.499886    0.500000   0.000114                    *
'* -----------------------------------------------                  *
'*                                                                  *
'*                              Basic Release By J-P Moreau, Paris. *
'*                                       (www.jpmoreau.fr)          *
'********************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

  H = .05  'integration step

  DIM B(4), X(4), Y(4)

  'Initial conditions
  X(0) = 0#      'starting X
  X1 = 1#        'ending X
  Y(0) = 1#      'initial Y 
  
  'write header
  CLS
  PRINT
  PRINT "       X           Y        Y exact      Error    "
  PRINT "  ------------------------------------------------"
  'write initial line
  F$ = "#####.######"
  PRINT USING F$; X(0);
  PRINT USING F$; Y(0);
  PRINT USING F$; Y(0);
  PRINT "    0"
  'use Runge-Kutta to start
  FOR K = 0 TO 1
    GOSUB 1000 'call RK4
    XX = X(K + 1): GOSUB 600: ER = ABS(FX - Y(K + 1))
    PRINT USING F$; X(K + 1);
    PRINT USING F$; Y(K + 1);
    PRINT USING F$; FX;
    PRINT USING F$; ER
  NEXT K
  'main integration loop
  WHILE X(2) < X1
    FOR I = 1 TO 3
      XX = X(3 - I): YY = Y(3 - I): GOSUB 500
      B(I) = F
    NEXT I
    X(3) = X(2) + H
    Y(3) = Y(2) + H * (23# * B(1) - 16# * B(2) + 5# * B(3)) / 12#
    XX = X(3): GOSUB 600: ER = ABS(Y(3) - FX)
    PRINT USING F$; X(3);
    PRINT USING F$; Y(3);
    PRINT USING F$; FX;
    PRINT USING F$; ER
    FOR K = 0 TO 2
      X(K) = X(K + 1): Y(K) = Y(K + 1)
    NEXT K
  WEND
  PRINT "  ------------------------------------------------"

END

'User defined function Y'=F(XX,YY)
500  F = -YY + XX / ((1# + XX) * (1# + XX))
     RETURN

'Exact solution Y=FX(XX)
600  FX = 1# / (1# + XX)
     RETURN

'Runge-Kutta method to calculate first points only
1000 XX = X(K): YY = Y(K): GOSUB 500: C1 = F
     XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GOSUB 500: C2 = F
     XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C2: GOSUB 500: C3 = F
     XX = X(K) + H: YY = Y(K) + H * C3: GOSUB 500: C4 = F
     X(K + 1) = X(K) + H
     Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6#
     RETURN

'end of file adambash.bas