'********************************************************************
'*             Differential equations of order N                    *
'*             by Runge-Kutta method of order 4                     *
'* ---------------------------------------------------------------- *
'* 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:                                                      *
'*                                                                  *
'* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6               *
'*          with initial conditions: y(4)=2 and y'(4)=-2tan(1)      *
'*                                                                  *
'* Exact solution is:   y = 2cos(1)/cos(x-5)                        *
'*                                                                  *
'*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N          *
'*              of type y(n) = f(y(n-1),...,y',y,x)                 *
'*                                                                  *
'*   order of equation: ? 2                                         *
'*   begin value x    : ? 4                                         *
'*   end value x      : ? 6                                         *
'*   y0 value at x0   : ? 2                                         *
'*   y1 value at x0   : ? -3.114815                                 *
'*   number of points : ? 11                                        *
'*   finesse          : ? 20                                        *
'*                                                                  *
'*       X            Y                                             *
'* ---------------------------                                      *
'*   4.000000     2.000000                                          *
'*   4.200000     1.551018                                          *
'*   4.400000     1.309291                                          *
'*   4.600000     1.173217                                          *
'*   4.800000     1.102583                                          *
'*   5.000000     1.080605                                          *
'*   5.200000     1.102583                                          *
'*   5.400000     1.173217                                          *
'*   5.600000     1.309291                                          *
'*   5.800000     1.551018                                          *
'*   6.000000     2.000000                                          *
'*                                                                  *
'********************************************************************
defint i-n
defdbl a-h,o-z

'ifi,i,ndata,iordre: INTEGER
dim yi(10), t(50)
dim ta(10),tb(10),tc(10),td(10),y(10),z(10)

  Cls
  print
  print "    DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N"
  print "          of type y(n) = f(y(n-1),...,y',y,x)"
  print
  print "    order of equation: "; : input iordre
  print
  print "    begin value x    : "; : input xi
  print "    end value x      : "; : input xf
  for i=0 to iordre-1
    print "    y";i;" value at x0   : "; : input yi(i)
  next i
  print "    number of points : "; : input m
  print "    finesse          : "; : input ifi

  'call subroutine equadiffn
  gosub 2000

END

'Example: y"=(2y'y'+yy)/y
1000 'FUNCTION fp(x,y())
  if abs(y(0))<1e-12 then y(0)=1e-12
  fp=(2*y(1)*y(1)+y(0)*y(0))/y(0)
return

'***************************************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER N        *
'*                of type y(n) = f(y(n-1),y(n-2),...,y',y,x)               *
'* ----------------------------------------------------------------------- *
'*  INPUTS:                                                                *
'*    fp        Equation to integrate                                      *
'*    xi, xf    Begin, end values of variable x                            *
'*    Yi        Begin values at xi (of f and derivatives)                  *
'*    m         Number of points to calculate                              *
'*    n         Order of differential equation                             *
'*    fi        finesse (number of intermediate points)                    *
'*                                                                         *
'*  OUTPUTS:                                                               *
'*    t         real vector storing m results for function y               *
'* ----------------------------------------------------------------------- *
'*  EXAMPLE:    y" = (2 y'y' + yy) / y with y(4) = 2, y'(4) = -2*tan(1)    *
'*              Exact solution:  y = (2 cos(1))/cos(x-5)                   *
'***************************************************************************
2000 'Subroutine Equadiffn
'h,x,a,b,c,d : double
'ta,tb,tc,td,y,z : Table;
'i,j,k,ni,n1,n2 : integer

   n=iordre
   if ifi<1 then return
   h = (xf - xi) / ifi / (m-1)
   n1=n-1 : n2=n-2
   t(1)=Yi(0)
   for k=0 to n1
     y(k)=Yi(k) : z(k)=Yi(k)
   next k
   for i=1 to m
     ni=(i-1)*ifi-1
     for j=1 to ifi
       x=xi+h*(ni+j)
       for k=0 to n1
         y(k)=z(k)
       next k
       gosub 1000
       a=h*fp
       for k=0 to n2
         ta(k)=h*y(k+1) : y(k)=z(k)+ta(k)/2#
       next k
       y(n1)=z(n1)+a/2#
       x=x+h/2
       gosub 1000
       b=h*fp
       for k=0 to n2
         tb(k)=h*y(k+1) : y(k)=z(k)+tb(k)/2#
       next k
       y(n1)=z(n1)+b/2#
       gosub 1000
       c=h*fp
       for k=0 to n2
         tc(k)=h*y(k+1) : y(k)=z(k)+tc(k)
       next k
       y(n1)=z(n1)+c
       x=x+h/2
       gosub 1000
       d=h*fp
       for k=0 to n2
         z(k)=z(k)+(ta(k)+2#*tb(k)+2#*tc(k)+h*y(k+1))/6#
       next k
       z(n1)=z(n1)+(a+b+b+c+c+d)/6#
     next j
     t(i+1)=z(0)
   next i
   'call subroutine Affiche
   gosub 3000
Return

3000 'Subroutine Affiche
  h=(xf-xi)/(m-1)
  x=xi-h
  Cls
  print
  print "      X         Y     "
  print "----------------------"
  for i=1 to m
    x=x+h
    print using " ##.######  ##.######"; x; t(i)
  next i
return

'End of file teqdifn.bas