'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
'* Reference:                                              *
'*   "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*    DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'*    [BIBLI 03].                                          *
'*                                                         *
'*                    Basic version by J-P Moreau, Paris   *
'*                            (www.jpmoreau.fr)            *
'* ------------------------------------------------------- *
'* SAMPLE RUN:                                             *
'* (Integrate y'=(4xy+4x*sqrt(y))/(1+xý) from x=0 to x=10  *
'*  with initial condition y(0)=1)                         *
'*                                                         *
'*  Input x begin: 0                                       *
'*  Input x end  : 10                                      *
'*  Input y begin: 1                                       *
'*  Input number of points: 10                             *
'*  Input finesse: 10                                      *
'*                                                         *
'*         X               Y                               *
'*  -----------------------------                          *
'*     1.000000         8.999984                           *
'*     2.000000        80.999881                           *
'*     3.000000       360.999496                           *
'*     4.000000      1088.998512                           *
'*     5.000000      2600.996484                           *
'*     6.000000      5328.992838                           *
'*     7.000000      9800.986875                           *
'*     8.000000     16640.977767                           *
'*     9.000000     26568.964559                           *
'*    10.000000     40400.946171                           *
'*  -----------------------------                          *
'*                                                         *                                                    *
'***********************************************************
defint i-n
defdbl a-h,o-z

  MAXDATA=100

  dim tx(MAXDATA),ty(MAXDATA)
  'integer ifi, n

  cls
  print
  input " Input x begin: ", xi
  input " Input x end  : ", xf
  input " Input y begin: ", yi
  input " Input number of points: ", n
  input " Input finesse: ", ifi

  gosub 3000  'call equadiff_pc(tx,ty,xi,xf,yi,n,ifi)

' print results
  F$="#####.######     #####.######"
  print
  print "       X                Y     "
  print " -----------------------------"
  for i=1 to n
    print using F$; tx(i); ty(i)
  next i
  print " -----------------------------"
  print

end 'of main program

1000  'Function y'=(4xy+4x*sqrt(y))/(1+xý)
  if yy>0 then
    f = (4*xx*yy+4*xx*sqr(yy))/(1+xx*xx)
  else
    f = 0
  end if
return

2000 'classical Runge-Kutta method of order 4
  xx=x : yy=y : gosub 1000 : a=h*f
  xx=x+h/2 : yy=y+a/2 : gosub 1000 : b=h*f
  yy=y+b/2 : gosub 1000 : c=h*f
  x=x+h : xx=x : yy=y+c : gosub 1000 : d=h*f
  y = y + (a+b+b+c+c+d)/6
return

  '**********************************************************
  '*             Prediction-correction method               *
  '* ------------------------------------------------------ *
  '* INPUTS:                                                *
  '*	       xi      begin x value                        *
  '*	       xf      end x value                          *
  '*           yi      begin y value (at x=xi)              *
  '*           n       number of points to calculate        *
  '*           ifi     finesse (number of intermediate      *
  '*                   points (for example 10)              *
  '* OUTPUTS:                                               *
  '*           tx      table of x values (m values)         *
  '*           ty      table of y values (m values)         *
  '*                                                        *
  '* DESCRIPTION:                                           *
  '* The algorithm has the following steps:                 *
  '*     1. calculate y2, y3, y4 using a classical Runge-   *
  '*        Kutta method of order 4.                        *
  '*     2. from point (x4, y4), first estimate y(n+1) by   *
  '*        formula:                                        *
  '*        y(n+1)=y(n) + h/24(55y'(n)-59y'(n-1)+37y'(n-2)  *
  '*               -9y'(n-3)                                *
  '*        then continue with formula:                     *
  '*        y(n+1)=y(n) + h/24(9y'(n+1)+19y'(n)-5y'(n-1)    *
  '*               +y'(n-2),                                *
  '*        noting that y'(n+1)=f(x(n+1),y(n+1)) with the   *
  '*        estimated value of y(n+1), until convergence is *
  '*        obtained.                                       *
  '**********************************************************
3000 'Subroutine equadiff_pc(tx, ty, xi, xf, yi, n, ifi)
  dim p(3)
  z=yi
  m=n
  if (m>MAXDATA) or (ifi<1) then return
  h = (xf-xi)/ifi/m
  xx=xi : yy=yi : gosub 1000 : p(3)=f
  tx(0)=xi : ty(0)=yi
  k=0
  for i=1 to m
    ni = (i-1)*ifi-1
    for j=1 to ifi
      x=xi+h*(ni+j)
      k=k+1
      if (k<4) then
        'call runge_kutta(h,x,z)
        y=z : gosub 2000 : z=y
	xx=x : yy=z : gosub 1000 : p(3-k)=f
      else
	x=x+h
        w=z+(h/24)*(55*p(0)-59*p(1)+37*p(2)-9*p(3))
        'continue
10	y=w
          xx=x : yy=y : gosub 1000
          w=z+(h/24)*(9*f+19*p(0)-5*p(1)+p(2))
        if abs(y-w) > 1e-10 then goto 10
	z=w : p(3)=p(2) : p(2)=p(1)
	p(1)=p(0) : xx=x : yy= z : gosub 1000 : p(0)=f
      end if
    next j 'j loop
    tx(i)=x : ty(i)=z
  next i 'i loop
return

' end of file teqdifpc.bas