  ╔══════════════════════════════╗                (C) Copyright 1986-1990
  ║  WELCOME to the ODE Solver ! ║                      Zvi Shippony
  ╚══════════════════════════════╝                     (818) 990-0134

    This Program solves Ordinary Differential Equations. The main use is
 solving Initial Condition (I.C) type problems but for second order equation
 the user have the option to solve a Two-Point-Boundary Value problem (B.V).

 Option #1:
 ----------
    Solves a set of N First order ODE whose Right-Hand-Side (RHS) is any
 legal expression of the form: F(X,Y0,Y1,Y2,...,Y9)   (Thus N <= 10 ).

 Option #2:
 ----------
    Solves an N-th order ODE whose RHS is any legal expression of the form:
 F(X,Y,Y',Y'',Y''',..,Y''''''''')  or F(X,Y0,Y1,Y2,Y3,..,Y9)  where Y0 stands
 for Y, Y1 stands for Y', Y2 stands for Y'' and so on. YOU CAN USE BOTH
 NOTATIONS.
$$$

 'Expression' is any legal combination of: +, -, *, /, **, (, ), AND any of
  the following functions:

   ABS, INT, EXP, SIN, COS, TAN, COT, LOG, LN, FACT or ! (Factorial), SIGN,
   SQRT, SINH, COSH, TANH, ARCSIN, ARCCOS, ARCTAN, ARCSINH, ARCCOSH, ARCTANH

   And the "Special Functions" :

     Z(x)     { Riemann's Zeta function }
     G(x)     { Gamma function, (IF x is an integer then x! = G(x+1)) }
     BJ(n,x)  { Bessel Function of the  first kind, J(n,x) }
     BY(n,x)  { Bessel Function of the second kind, Y(n,x) }
     BI(n,x)  { Modified Bessel Function of the  first kind, I(n,x) }
     BK(n,x)  { Modified Bessel Function of the second kind, K(n,x) }
     SBJ(n,x) { Spherical Bessel Function of the  first kind, j(n,x) }
     SBY(n,x) { Spherical Bessel Function of the second kind, y(n,x) }

  PI is a reserved name and will be interpeted as Pi = 3.14159265358....

$$$
    Interval of solution (A,B) is needed, where A  and/or B  are either given
 as numbers or any expression containing numbers. A comma is required between
 A and B .

    For example: if A = 0.0 and B = Pi/4.0, the input could be given as:
       0.0 , pi/4   or:  0.0 , arctan(1.0)   .

    Initial Conditions (or Boundary Values) could be given as numbers or any
 legal expression containing numbers .

 NOTE ON TWO-POINT-BOUNDARY-VALUE PROBLEMS: (Applies for 2-nd order eqn. only)
 -----------------------------------------------------------------------------
    Only explicit boundary values are accepted. That is, boundary  value can
 not be a function of Y itself.
   ( Right example: Y(1.0)=5*exp(1)    Wrong example: Y(1.0)=5*Y(0.0)  )

    The method for solving boundary value problem is a combination of FINITE-
 DIFFERENCE Method and the SHOOTING Method.  The user can supply the program
 with the partial derivatives of F(X,Y,Y') w.r.t.  Y and Y', this way the
 convergence will be faster (using Newton's Method) but if these are not
 supplied, the iterations will be done with the False-Position Method, which
 gives slower convergence (and the FINITE-DIFFERENCE step will NOT be used).
$$$
  Accuracy of solution:
  ---------------------
    The program will try to achieve the inputted accuracy, but the higher
 the requested accurate is,  the longer the execution time will be. Some time
 roundoff errors will stop the execution without achieving the desired
 accuracy. The smallest possible accuracy is: 5.0E-14 (default).

    The method used for all I.C. problems is the Bulirsch_Stoer integration
 procedure with 4-th order Runge-Kutta step taken at points of non-convergence.

    When the program writes: '** Finished' , the solution is stored in an
 array of up to 100 entries equally spaced.

    The user has the option to ask for any specific X-value (as long as this
 value is in the interval [A,B] ) and the program will interpolate for that
 value.

    Another way is to store these entries in a file for later use (Plotting,
 Printing,  etc.)
$$$

 Option #3:  Loading a system from a file
 ----------
    Any ASCII file written in the format bellow could be used to load a pre-
 determined system into the program for solving. The file MUST have that
 format in order to work properly. It could be written in any editing type
 S/W (like: Edlin, WordStar, Word-Perfect, Word  etc..). Each line has a colon
 (':') separating between the name of the input and its value.

 Note: You can insert comments anywhere in the line as { ... }.


 There are 2 slightly different file formats, according to the option needed:
$$$

         File format for option 1 (System of N first order equations):
         -------------------------------------------------------------

Option  : 1
Order   : {1 to 10}
Type    : IC
Domain  : {X-Range in the format: x1 , x2}  (Don't forget the comma !!)
Accuracy: {Integration accuracy to use}
StepSize: {Step size to use}
Eqn. #0 : {The Right-hand-size expression to use for Y0'}
Eqn. #1 : {The Right-hand-size expression to use for Y1'}
Eqn. #2 : {The Right-hand-size expression to use for Y2'}
   .                    .
   .                    .
Eqn. #n :     {Where n is = order -1 }
Initial condition Y0(LHS): {Y  value at the initial point}
Initial condition Y1(LHS): {Y' value at the initial point}
   .                                    .
   .                                    .
Initial condition Yn(LHS):     {Where n is = order-1}

$$$
                File format for option 2 (N-th order equation):
                -----------------------------------------------

Option  : 2
Order   : {1 to 10}
Type    : {IC or BV}
Domain  : {X-Range in the format: x1 , x2}  (Don't forget the comma !!)
Accuracy: {Integration accuracy to use}
StepSize: {Step size to use}
Equation: {The Right-hand-size expression to use}
dF/dY   : {FOR BV PROBLEMS ONLY, the derivative of F(X,Y,Y') w.r.t  Y }
dF/dY'  : {FOR BV PROBLEMS ONLY, the derivative of F(X,Y,Y') w.r.t  Y'}
Boundary Value Y0(LHS) : {The Left-Hand_Side Y value }
Boundary Value Y0(RHS) : {The Right-Hand-Side Y value}
       or:
Initial condition Y0(LHS): {Y  value at the initial point}
Initial condition Y1(LHS): {Y' value at the initial point}
          .                              .
          .                              .
Initial condition Yn(LHS):     {Where n is the order-1}
$$$

NOTE:

    Some times the equation (the 7-th input) is too long to fit on one line.
 In this case you can continue it on the next line by starting it with
 a colon (':') and then type in the rest of the equation, i.e. :

Equation: (1 + sin(y') - cos(y''*y''') + tan(y''''*y''''') -
 : y''*y'''') / sqrt(y''**2 + y'''**2)


  YOU CAN USE AS MANY CONTINUATIONS YOU WANT,BUT AFTER SQUEEZING ALL THE
  EXPRESSION(s) FROM BLANKS AND CONVERTING ALL THE PRIME NOTATIONS
  (Y',Y'',Y''',...) INTO (Y1,Y2,Y3,..) THE FINAL EQUATION CAN NOT BE LONGER
  THEN 76 CHARACTERS !!

$$$

 Option #4:  Storing a system into a file
 ----------
    After a system has been solved, you can store it into a file. The program
 will write an ASCII file in the format given in the above example. This file
 can be then edited and used later as loading file (as in option 3 above)

$$$




 WE ARE GOING TO WALK THROUGH 2 COMPLETE EXAMPLES NOW ..

$$$

  Example #1:    After the program is loaded, it starts like:
  -----------

  Your options are:

  1 - Solve system of N first order O.D.E.: Yi' = F(X,Yj) (i,j = 0,1,..,N-1)
  2 - Solve N-th order O.D.E., e.g: Y''' = F(X,Y0,Y1,Y2)
      Where: Y0 ≡ Y, Y1 ≡ Y', Y2 ≡ Y''   (You can use either notation..)
  3 - Load a system from a file
  4 - Store last system on file
  M - Memory manipulations
  B - Set Bell  ON / OFF  for error or info. msg. Bell is now: OFF
  H - Help and information about the program
    - Any other key to quit ..

  Enter your Choice (1 - 4, M, B, H or Quit):

 ( LET US CHOOSE 2nd-ORDER  2-POINTS BOUNDARY-VALUE PROBLEM: Y''= -(2*Y'+Y)
   Solution: Y(x)=X*EXP(-X) , FOR X=[-1,4], Y(-1)=-EXP(1), Y(4)=4*EXP(-4) )

  Enter your Choice (1, 2, 3, B, H or Quit): 2
  Store various inputs in Memory ? (y/n): n   (SAY: N)
  Enter order of ODE (1-10): 2
$$$
   Enter: IC for Initial Condition Problem
   Enter: BV for Two-Point-Boundary Value Problem
          Your choice: BV

 (THE PROGRAM WILL NOW PRESENTS YOU WITH THE LINE: Y'' =
  THIS IS THE TIME TO ENTER THE EQUATION(S) ...)

  Y'' = -(2*y'+y)        ( SO F(X,Y,Y') = -2*Y'-Y )
 (YOU COULD ALSO ENTER: -(2*Y1+Y0)  or: -(2*Y'+Y0) etc ..)

** For faster convergence we need dF/dY and dF/dY'
   Are you supplying those ? (y/n): y               (SAY: Y)

  dF/dY  = -1                   (THIS IS DF/DY )
  dF/dY' = -2                   (THIS IS DF/DY')

  Enter x Domain ( e.g.:  -1.0 , 1.0): -1,4
  Enter desired accuracy (Default is: 5.0E-14): 1.e-10
  Compute solution every what Delta-x ? 0.05           (Total of 100 points)

  Enter B.V:  Y(-1.0) (Default: 0.0): -exp(1)
  Enter B.V:  Y(4.0) (Default: 0.0): 4*exp(-4)
  ** Do you wish to stop (pause) after each iteration ? (y/n): n  (SAY: N )
 ( YOU HAVE THE OPTION TO CONTROL THE ITERATIONS PROCESS )
$$$
 ( SINCE WE CHOSE 2-POINT BOUNDARY-VALUE PROBLEM, THE FIRST METHOD IS THE
   FINIT-DIFFERENCE METHOD ..THE PROGRAM WILL SAY:)

** Iterative Finite-Difference Method  **
   # of iterations:   2 Yields: dy/dx(LHS) =  5.43638706110044E+000
   Number of mesh points: 100  Step size: 4.95050E-002
   Solution accuracy at mesh points: Order(2.0E-3)

   Are you satisfied with the current accuracy ? (y/n) n     (SAY: N )

   Start Shooting Method now with the above guess for dy/dx(LHS)

                                ( Working ... )
$$$

** Iterative Shooting Method  **
   Given B.V:  7.32625555549367E-002  Achieved B.V:  7.32625555469740E-002
   Rel. Error: 7.96270E-012   Using dy/dx(LHS) =  5.43656365665446E+000
   Integration accuracy: 1.00000E-010  Max. step size used: 5.00000E-002
   Next guess for dy/dx(LHS)  is:  5.43656365668701E+000

                                ( Working ... )

  *** Finished ...        (HOORAY ! IT CONVERGED ...)

  To Plot or do Least-Square-Fit of the solution we need to write the results
  into a file.  Write a file ? (y/n):    (SAY: Y)
  Enter file name: ODETST

First title for the file is: SOLUTION OF: Y'' = -(2*Y1+Y0)
Enter (optional) 2-nd title (Initial or Boundary conditions, remarks, etc ..)
Title:Testing. BV: y(-1)=-exp(1), y(4)=4*exp(-4)  eps=1.0e-10
(THE PROGRAM WILL RESPOND BY:) ** File: ODETST is ready .
$$$
  Example #2:
  -----------
  Your options are:

  1 - Solve system of N first order O.D.E.: Yi' = F(X,Yj) (i,j = 0,1,..,N-1)
  2 - Solve N-th order O.D.E., e.g: Y''' = F(X,Y0,Y1,Y2)
      Where: Y0 ≡ Y, Y1 ≡ Y', Y2 ≡ Y''   (You can use either notation..)
  3 - Load a system from a file
  4 - Store last system on file
  M - Memory manipulations
  B - Set Bell  ON / OFF  for error or info. msg. Bell is now: OFF
  H - Help and information about the program
    - Any other key to quit ..

  Enter your Choice (1 - 4, M, B, H or Quit):
 ( LET US CHOOSE SYSTEM OF TWO 1-st ORDER EQUATIONS)

  Enter your Choice (1, 2, 3, B, H or Quit): 1
  Store various inputs in Memory ? (y/n): n   (SAY: N)
  Enter number of equations (1-10): 2

  Y0' = Y1
  Y1' = -Y0      ( The solution is: Y0(X)=SIN(X), Y1(X)=COS(X)  )
$$$

  Enter x Domain ( e.g.:  -1.0 , 1.0): -PI,PI
  Enter desired accuracy (Default is: 5.0E-14): 1.E-8
  Compute solution every what Delta-x ? 0.1

  Enter I.C:  Y0(-3.14159265358979) (Default: 0.0): 0.0
  Enter I.C:  Y1(-3.14159265358979) (Default: 0.0): -1.0

                                ( Working ... )

  *** Finished ...

  To Plot or do Least-Square-Fit of the solution we need to write the results
  into a file.  Write a file ? (y/n):    (SAY: N)
       ( WE ARE SKIPPING THE FILE WRITING OPTION .. )
$$$
( THE PROGRAM WILL SAY:)

  For what x you want y(x) ? (Hit Enter to stop) 1.1    (ENTER: 1.1 )

( THE PROGRAM WILL INTERPOLATE THE SOLUTION FOR X = 1.1 )

 x:   1.10000000000000E+000 y0:   8.91207360381939E-001
                            y1:   4.53596121602802E-001

( CHECK IT OUT, THE TRUE SOLUTION IS:

   SIN(1.1) = 8.91207360061435E-001   SO THE ERROR IS:  3.21e-11
   COS(1.1) = 4.53596121425577E-001   SO THE ERROR IS:  1.77e-11  )

( THE PROGRAM WILL CONTINUE TO ASK FOR WHAT X TO INTERPOLATE, HIT ENTER TO
  QUIT ..)

( THE PROGRAM WILL GIVE YOU ONE MORE CHANCE TO WRITE THE SOLUTION INTO A FILE
  AS DESCRIBED IN EXAMPLE #1 ABOVE ..)

                            That's all folks ...
