'********************************************************************
'*             Differential equations of order 1                    *
'*             by Runge-Kutta method of order 4                     *
'*                                                                  *
'*                             Basic version by J-P Moreau, Paris   *
'*                                     (www.jpmoreau.fr)            *
'* ---------------------------------------------------------------- *
'* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'*             [BIBLI 03].                                          *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'*                                                                  *
'* Example: integrate y'=4x*(y+sqrt(y))/1+x^2 from x=0 to x=1       *
'*                                                                  *
'*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1          *
'*                     of type y' = f(x,y)                          *
'*                                                                  *
'*   begin value x   : ? 0                                          *
'*   end value x     : ? 1                                          *
'*   y value at x0   : ? 1                                          *
'*   number of points: ? 11                                         *
'*   finesse         : ? 10                                         *
'*                                                                  *
'*       X            Y                                             *
'* -------------------------                                        *
'*   0.000000     1.000000                                          *
'*   0.100000     1.040400                                          *
'*   0.200000     1.166400                                          *
'*   0.300000     1.392400                                          *
'*   0.400000     1.742400                                          *
'*   0.500000     2.250000                                          *
'*   0.600000     2.958400                                          *
'*   0.700000     3.920400                                          *
'*   0.800000     5.198400                                          *
'*   0.900000     6.864400                                          *
'*   1.000000     9.000000                                          *
'*                                                                  *
'********************************************************************
defint i-n
defdbl a-h,o-z

'ifi:INTEGER

  Cls
  print
  print "    DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1"
  print "             of type y' = f(x,y)"
  print
  print "    begin value x   : "; : input xi
  print "    end value x     : "; : input xf
  print "    y value at x0   : "; : input yi
  print "    number of points: "; : input m
  print "    finesse         : "; : input ifi

  dim t(m+1)

  'equadiff1(fp,t,xi,xf,yi,m,ifi)
  gosub 2000

END


1000 'Example: y'=4x(y+rac(y))/(1+x²)
    fp=4#*xx*(yy+sqr(yy))/(1#+xx*xx)
return

'***************************************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER 1        *
'*                       of type y' = f(x,y)                               *
'* ----------------------------------------------------------------------- *
'*  INPUTS:                                                                *
'*    fp        Given equation to integrate (see test program)             *
'*    xi, xf    Begin, end values of variable x                            *
'*    yi        Begin Value of y at x=xi                                   *
'*    m         Number of points to calculate                              *
'*    fi        finesse (number of intermediate points)                    *
'*                                                                         *
'*  OUTPUTS:                                                               *
'*    t       : real vector storing m results for function y               *
'***************************************************************************
2000 'Subroutine  Equadiff1
  if ifi < 1 then return
  h   = (xf - xi) / ifi / (m-1)
  y   = yi
  t(1) = yi

  for i = 1 to m
    ni = (i - 1) * ifi - 1
    for j = 1 to ifi
      x = xi + h * (ni + j)
      xx=x : yy=y
      gosub 1000 'calculate current fp(x,y)
      a = h * fp
      xx=x+h/2# : yy=y+a/2#
      gosub 1000 'calculate fp(x+h/2,y+a/2)
      b = h * fp
      yy=y+b/2#
      gosub 1000 'calculate fp(x+h/2,y+b/2)
      c = h * fp
      x = x + h
      xx=x : yy=y+c
      gosub 1000 'calculate fp(x,y+c)
      d = h * fp
      y = y + (a + b + b + c + c + d) / 6
    next j
    t(i+1) = y
  next i
  'Affiche(t,m,xi,xf)
  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 teqdif1.bas