!************************************************************************
!*                                                                      *
!* Test examples for methods to solve first order ordinary systems of   *
!* differential equations.                                              *
!*                                                                      *
!* We supply the right hand sides, their alpha-numeric description and  *
!* the exact analytic solution, if known.                               *
!*                                                                      *
!* When running the main test program, the user can select among the    *
!* examples below.                                                      *
!*                                                                      *
!*                               F90 1.1 Release By J-P Moreau, Paris.  *
!* -------------------------------------------------------------------- *
!* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges        *
!*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
!*                                                                      *
!* Release 1.1: added example #1 for test program tawp.                 *
!************************************************************************

! Note: here, only examples #0, #1 and #3 are implemented in Fortran
  
  Subroutine dgl(bspn,n,x,y,f)
    real*8 x,y(0:n-1), f(0:n-1)
    integer bspn, n
    Select Case(bspn)
      Case (0)
        f(0) = y(0) * y(1) + DCOS(x) - 0.5D0 * DSIN(2.d0 * x)
        f(1) = y(0)**2 + y(1)**2 - (1.d0 + DSIN(x))
      Case (1)
	f(0) = -y(0) + x / ((1.d0+x)*(1.d0+x))
      Case (3)
        f(0) = y(1)
        f(1) = y(2)
        f(2) = y(3)
        f(3) = y(4)
        f(4) = (45.d0 * y(2) * y(3) * y(4) - 40.d0 * y(3)**3) / (9.d0 * y(2)**2)
    End Select
  return
  end

  Subroutine settxt(bspn,n,dgltxt)
    integer bspn, n
	character*70 dgltxt(0:n-1)
    Select Case(bspn)
      Case (0)
        n=2
        dgltxt(0) = ' y1'' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)'
        dgltxt(1) = ' y2'' = y1 * y1 + y2 * y2 - (1 + sin(x))'
      Case (1)
        n=1
	dgltxt(0) = ' y'' = -y + x / ((1+x)*(1+x))'
      Case (3)
        n=5
        dgltxt(0) = ' y1'' = y2'
        dgltxt(1) = ' y2'' = y3'
        dgltxt(2) = ' y3'' = y4'
        dgltxt(3) = ' y4'' = y5'
        dgltxt(4) = ' y5'' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)'
    End Select
  return
  end

!Other examples given in C

! ----------------------- DE system number 0 -----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl0(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(0) * y(1) + COS(x) - HALf * SIN(TWO * x);
!  f(1) = y(0) * y(0) + y(1) * y(1) - (ONE + SIN(x));
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt0(void)
!{
!  return
!    "y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)\n"
!    "y2' = y1 * y1 + y2 * y2 - (1 + sin(x))\n";
!}


! ----------------------- DE system number 1 -----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl1(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(1);
!  f(1) = (ONE + y(1) * y(1)) / y(0);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt1(void)
!{
!  return
!    "y1' = y2\n"
!    "y2' = (1 + y2 * y2) / y1\n";
!}


! ----------------------- DE system number 2 -----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl2(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(1);
!  f(1) = -fOUR * y(0) + EXP(x);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt2(void)
!{
!  return
!    "y1' = y2\n"
!    "y2' = -4 * y1 + exp(x)\n";
!}


!----------------------- DE system number 3 ------------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl3(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(1);
!  f(1) = y(2);
!  f(2) = y(3);
!  f(3) = y(4);
!  f(4) = ((REAL)45.0 * y(2) * y(3) * y(4) -
!          (REAL)40.0 * y(3) * y(3) * y(3)) / (NINE * y(2) * y(2));
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt3(void)
!{
!  return
!    "y1' = y2\n"
!    "y2' = y3\n"
!    "y3' = y4\n"
!    "y4' = y5\n"
!    "y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)\n";
!}



!----------------------- DE system number 4 ------------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl4(REAL x, REAL *y, REAL *f)
!{
!  f(0) = -y(1);
!  f(1) = y(0);
!  f(2) = y(4) - y(1) * y(2);
!  f(3) = y(0) * y(3) - y(5);
!  f(4) = -y(2) - y(1) * y(4);
!  f(5) = y(3) + y(0) * y(5);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt4(void)
!{
!  return
!    "y1' = -y2\n"
!    "y2' = y1\n"
!    "y3' = y5 - y2 * y3\n"
!    "y4' = y1 * y4 - y6\n"
!    "y5' = -y3 - y2 * y5\n"
!    "y6' = y4 + y1 * y6\n";
!}

!*************************************************************************
!* Compute the value of the analytic solution y(x) for the above DE      *
!* for the initial value problem  y(0) = (1,0,0,1,e,0)                   *
!*************************************************************************
!static void loesung4(REAL x, REAL *y)
!{
!  y(0) = COS(x);
!  y(1) = SIN(x);
!  y(2) = SIN(x) * EXP(COS(x));
!  y(3) = COS(x) * EXP(SIN(x));
!  y(4) = COS(x) * EXP(COS(x));
!  y(5) = SIN(x) * EXP(SIN(x));
!}


!--------------------------- DE system number 5 --------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl5(REAL x, REAL *y, REAL *f)
!{
!  f(0) = (REAL)-500.5 * y(0) + (REAL)499.5 * y(1);
!  f(1) =  (REAL)499.5 * y(0) - (REAL)500.5 * y(1);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt5(void)
!{
!  return
!    "y1' = -500.5 * y1 + 499.5 * y2\n"
!    "y2' =  499.5 * y1 - 500.5 * y2\n";
!}

!*************************************************************************
!* Compute the value of the analytic solution y(x) for the above DE      *
!* for the initial value problem  y(0) = (4,2)                           *
!*************************************************************************
!static void loesung5(REAL x, REAL *y)
!{
!  y(0) = THREE * EXP(-x) + EXP((REAL)-1000.0 * x);
!  y(1) = THREE * EXP(-x) - EXP((REAL)-1000.0 * x);
!}


!------------------------- DE system number 6 ----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl6(REAL x, REAL *y, REAL *f)
!{
!  f(0) = (REAL)15.0 * y(1);
!  f(1) = (REAL)-0.6 * y(0);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt6(void)
!{
!  return
!    "y1' = 15.0 * y2\n"
!    "y2' = -0.6 * y1\n";
!}

!*************************************************************************
!* Compute the value of the analytic solution y(x) for the above DE      *
!* for the initial value problem  y(0) = (0,1)                           *
!*************************************************************************
!static void loesung6(REAL x, REAL *y)
!{
!  y(0) = fIVE * SIN(THREE * x);
!  y(1) = COS(THREE * x);
!}



!------------------------- DE system number 7 ----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl7(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(1);
!  f(1) = y(2);
!  f(2) = y(3);
!  f(3) = fOUR * y(3) - SIX * y(2) + fOUR * y(1) - y(0);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt7(void)
!{
!  return
!    "y1' = y2\n"
!    "y2' = y3\n"
!    "y3' = y4\n"
!    "y4' = 4.0 * y4 - 6.0 * y3 + 4.0 * y2 - y1\n";
!}

!*************************************************************************
!* Compute the value of the analytic solution y(x) for the above DE      *
!* for the initial value problem  y(0) = (0,0,0,6)                       *
!*************************************************************************
!static void loesung7(REAL x, REAL *y)
!{
!  y(0) = x * x * x * EXP(x);
!  y(1) = x * x * (x + THREE) * EXP(x);
!  y(2) = x * (x * x + SIX * x + SIX) * EXP(x);
!  y(3) = (x * x * x + NINE * x * x + (REAL)18.0 * x + SIX) * EXP(x);
!}


!------------------------- DE system number 8 ----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl8(REAL x, REAL *y, REAL *f)
!{
!  f(0) = (REAL)-52.0 * y(0) + (REAL)50.0 * y(1);
!  f(1) =  (REAL)50.0 * y(0) - (REAL)52.0 * y(1);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt8(void)
!{
!  return
!    "y1' = -52 * y1 + 50 * y2\n"
!    "y2' =  50 * y1 - 52 * y2\n";
!}

!*************************************************************************
!* Compute the value of the analytic solution y(x) for the above DE      *
!* for the initial value problem                                         *
!*      y(0) = (1265588.55375228,-1265586.55375225)                      *
!*************************************************************************
!static void loesung8(REAL x, REAL *y)
!{
!  y(0) = (REAL)1.000000015       * EXP(        -TWO * x) +
!         (REAL)1265587.553752265 * EXP((REAL)-102.0 * x);
!  y(1) = (REAL)1.000000015 *       EXP(        -TWO * x) -
!         (REAL)1265587.553752265 * EXP((REAL)-102.0 * x);
!}


!------------------------- DE system number 9 ----------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl9(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(1);
!#if defined(BC3)
!  _fpreset();
!#endif
!  f(1) = -y(0) * COSH(x);
!}


!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt9(void)
!{
!  return
!    "y1' = y2\n"
!    "y2' = -y1 * cosh(x)\n";
!}


!------------------------- DE system number 10  --------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl10(REAL x, REAL *y, REAL *f)
!{
!  f(0) = y(1);
!  f(1) = -y(0) * y(0) * y(0);
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt10(void)
!{
!  return
!    "y1' = y2\n"
!    "y2' = -y1^3\n";
!}


!------------------------- DE system number 11 ---------------------------

!*************************************************************************
!* Compute the value f of the right hand side of a DE system at (x,y)    *
!*************************************************************************
!static void dgl11(REAL x, REAL *y, REAL *f)
!{
!#define dk  ((REAL)-10000.0)
!  f(0) = y(1);
!  f(1) = dk * y(0) + (dk - ONE) * y(1);
!#undef  dk
!}

!*************************************************************************
!* alpha-numeric description of above system                             *
!*************************************************************************
!static char *dgltxt11(void)
!{
!  return
!    "y1' =  y1                                    (dk = -10000)\n"
!    "y2' =  dk * y1 + (dk - 1) * y2\n";
!}

!*************************************************************************
!* Compute the value of the analytic solution y(x) for the above DE      *
!* for the initial value problem  (y0) = (1,1)                           *
!*************************************************************************
!static void loesung11(REAL x, REAL *y)
!{
!  y(0) = EXP(-x);
!  y(1) = -EXP(-x);
!}

!end of file t_dgls.f90