      program main

c*********************************************************************72
c
cc mus_test() tests mus().
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solution of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
      call timestamp ( )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'mus_test():'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test mus().'

      call musl_test ( )
      call musn_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'mus_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ''
      call timestamp ( )

      stop
      end
      subroutine musl_test ( )

c*********************************************************************72
c
cc musl_test() tests a linear problem.
c
c  Discussion:
c
c    Consider the ordinary differential equation of the general form:
c
c      dx(t) / dt = L(t)*x(t) + r(t),   0 <= t <= 6
c
c    and a boundary condition of the general form:
c
c      MA x(0) + MB x(6) = BCV  
c
c    where, for this problem, we have:
c
c      0 <= t <= 6
c
c             | 1-2cos(2t)    0   1+2sin(2t) |
c      L(t) = |   0           2     0        |
c             | -1+2sin(2t)   0   1+2cos(2t) |
c
c             | (-1+2cos(2t)-2sin(2t)) exp(t) |
c      r(t) = |       -exp(t)                 |
c             | ( 1-2cos(2t)-2sin(2t)) exp(t) |
c
c      MA = MB = I
c
c            | 1+exp(6) |
c      BCV = | 1+exp(6) |
c            | 1+exp(6) |
c 
c    The solution of this problem is:
c 
c             | exp(t) |
c      x(t) = | exp(t) |
c             | exp(t) |
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solution of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
      implicit none

      integer n
      parameter ( n = 3 )
      integer nti
      parameter ( nti = 15 )
      integer nu
      parameter ( nu = 6 )

      integer liw
      parameter ( liw = 3 * n )
      integer lw
      parameter ( lw = 8 * n + 2 * n * n )

      double precision a
      double precision ae
      double precision amp
      double precision b
      double precision bcv(n)
      double precision d(n,nti)
      double precision er(5)
      double precision exsol
      external fdifl
      external flin
      integer i
      integer ierror
      integer ihom
      integer iw(liw)
      integer j
      integer k
      integer kpart
      double precision ma(n,n)
      double precision mb(n,n)
      integer nrti
      double precision phirec(nu,nti)
      double precision q(n,n,nti)
      double precision ti(nti)
      double precision u(nu,nti)
      double precision w(lw)
      double precision x(n,nti)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'musl_test()'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test a linear problem.'
c
c  Set input parameters.
c
      ihom = 1
      ER(1) = 1.0D-11
      ER(2) = 1.0D-06
      ER(3) = 1.1D-16
      a = 0.D+00
      b = 6.D+00
      amp = 0.0D+00
      nrti = 10
c
c  Set the boundary condition matrices MA and MB.
c
      DO I = 1 , N
        DO J = 1 , N
          MA(I,J) = 0.0D+00
          MB(I,J) = 0.0D+00
        end do
        MA(I,I) = 1.D+00
        MB(I,I) = 1.D+00
      end do
c
c  Set the boundary condition vector BCV.
c
      BCV(1) = 1.0D+00 + DEXP(6.0D+00)
      BCV(2) = BCV(1)
      BCV(3) = BCV(1)
c
c  Call MUSL
c
      CALL MUSL ( FLIN, fdifl, N, IHOM, A, B, MA, MB, BCV, AMP, ER, 
     &  NRTI, TI, NTI, X, U, NU, Q, D, KPART, PHIREC, W, LW, IW, LIW, 
     &  IERROR )

      IF ( ( IERROR .NE. 0 ) .AND.
     &     ( IERROR .NE. 200 ) .AND. 
     &     ( IERROR .NE. 213 ) .AND.
     &     ( IERROR .NE. 240 ) ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'musl_test(): Fatal error!'
        write ( *, '(a,i6)' ) '  Error code IERROR = ', ierror
        write ( *, '(a)' ) '  Abnormal end of execution.'
        write ( *, '(a)' ) ''
        call timestamp ( )
        stop
      end if
c
c  Compare the computed and exact solutions.
c
      write ( *, '(a)' ) ''
      write ( *, '(a,g14.6)' ) '  Condition number =     ', er(4)
      write ( *, '(a,g14.6)' ) '  Amplification factor = ', er(5)
      write ( *, '(a)' ) ' '
      WRITE ( *, '(a)' ) '  I       T        APPROX. SOL.         '
     &  // 'EXACT SOL.        ABS. ERROR'
      write ( *, '(a)' ) ''
      DO K = 1 , NRTI
        EXSOL = DEXP(TI(K))
        AE = EXSOL - X(1,K)
        WRITE(*,'(i3,3x,f7.4,3(3x,d16.9))') K,TI(K),X(1,K),EXSOL,AE
        DO I = 2 , N
          AE = EXSOL - X(I,K)
          WRITE(*,'(13x,3(3x,d16.9))') X(I,K),EXSOL,AE
        end do
      end do
c
c  Terminate.
c
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'musl_test():'
      write ( *, '(a)' ) '  Normal end of execution.'

      return
      end
      subroutine flin ( t, y, f )

c*********************************************************************72
c
cc FLIN evaluates the linear part of the ODE, the value of L(t)*y(t).
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solutioni of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
c  Parameters:
c
c    Input, double precision T, the location in [A,B].
c
c    Input, double precision Y(N), the solution value at T.
c
c    Output, double precision F(N), the value of L(T) * Y(T).
c
      implicit none

      double precision co
      double precision f(3)
      double precision si
      double precision t
      double precision ti
      double precision y(3)

      TI = 2.0D+00 * T
      SI = 2.0D+00 * DSIN(TI)
      CO = 2.0D+00 * DCOS(TI)

      F(1) = (1.0D+00 - CO) * Y(1) + (1.0D+00 + SI) * Y(3)
      F(2) = 2.0D+00 * Y(2)
      F(3) = (-1.0D+00 + SI) * Y(1) + (1.0D+00 + CO) * Y(3)

      return
      end
      subroutine fdifl ( t, y, f )

c*********************************************************************72
c
cc fdifl evaluates the right hand side of the ODE, L(t)*y(t)+r(t).
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solutioni of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
c  Parameters:
c
c    Input, double precision T, the location in [A,B].
c
c    Input, double precision Y(N), the solution value at T.
c
c    Output, double precision F(N), the value of L(T) * Y(T) + R(T).
c
      implicit none

      double precision co
      double precision f(3)
      double precision si
      double precision t
      double precision ti
      double precision y(3)
c
c  Save time and have FLIN compute the linear part.
c
      CALL FLIN ( T, Y, F )
c
c  Add on the nonlinear part.
c
      TI = 2.0D+00 * T
      SI = 2.0D+00 * DSIN ( TI )
      CO = 2.0D+00 * DCOS ( TI )
      TI = DEXP ( T )

      F(1) = F(1) + ( -1.0D+00 + CO - SI ) * TI
      F(2) = F(2) - TI
      F(3) = F(3) + ( 1.0D+00 - CO - SI ) * TI

      RETURN
      END
      subroutine musn_test ( )

c*********************************************************************72
c
cc musn_test() tests a nonlinear problem.
c
c  Discussion:
c
c    The problem to be solved is as follows:
c
C          X1' = A * X1 / X2 * (X3 - X1 )
C          X2' = - A * ( X3 - X1 )
C          X3' = 1 / X4 * ( B - C * ( X3 - X5 ) - A * X3 * ( X3 - X1 ) )
C          X4' = A * ( X3 - X1 )
C          X5' = - C / D * ( X5 - X3 )
C
C          X1(0) = 1
C          X2(0) = 1
C          X3(0) = 1
C          X4(0) = -10
C          X3(1) = X5(1)
C
C          A = 0.5
C          B = 0.9
C          C = 1000.0
C          D = 10.0
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solution of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
c    Michael Hanke, Rene Lamour, Renate Winkler,
c    The program system RWA for the solution of two-point 
c    boundary-value problems: foundations, algorithms, comparisons,
c    Seminarberichte 67,
c    Sektion Mathematik, Humboldt University.
c
      implicit none

      integer lwg
      parameter ( lwg = 20 )
      integer n
      parameter ( n = 5 )
      integer nti
      parameter ( nti = 12 )
      integer nu
      parameter ( nu = 15 )

      integer liw
      parameter ( liw = 3 * n + nti )
      integer lw
      parameter ( lw = 7 * n + 3 * n * nti + 4 * n * n )

      double precision a
      double precision amp
      double precision b
      double precision d(n,nti)
      double precision er(5)
      external fdifn
      external g
      integer ierror
      integer itlim
      integer iw(liw)
      integer j
      integer k
      integer kpart
      integer nrti
      double precision phirec(nu,nti)
      double precision q(n,n,nti)
      double precision ti(nti)
      double precision u(nu,nti)
      double precision w(lw)
      double precision wgr(lwg)
      double precision x(n,nti)
      external x0t

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'musn_test():'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test a nonlinear problem.'
c
c  Set input parameters.
c
      ER(3) = 1.1D-15
      ER(1) = 1.0d-06
      ER(2) = 1.0D-02
      A = 0.0D+00
      B = 1.0D+00
      NRTI=10
      AMP=100.0D+00
      ITLIM=20
      IERROR=1
c
c  Call MUSN for the solution.
c
      CALL MUSN(fdifn,X0T,G,N,A,B,ER,TI,NTI,NRTI,AMP,ITLIM,X,Q,U,NU,D,
     1          PHIREC,KPART,W,LW,IW,LIW,WGR,LWG,IERROR)

      WRITE(*,*) ' MUSN: IERROR =',IERROR

      WRITE(*,200) A,B,ER(1),ER(2),ER(4),ER(5),KPART
  200 FORMAT(' A = ',F8.4,3X,'B = ',F8.4,/,' REQUIRED TOLERANCE = ',1P,
     1 D12.5,3X,'START TOLERANCE  = ',D12.5,/,
     2 ' CONDITION NUMBER = ',D12.5,3X,
     3 'AMPLIFICATION FACTOR = ',D12.5,/,' K-PARTITIONING =',I2,/)

      IF (IERROR.NE.0) GOTO 3000

      WRITE(*,215)
  215 FORMAT(' I ',4X,'T',9X,'Y1',12X,'Y2',12X,'Y3',12X,'Y4',12X,'Y5',
     1 /)
      DO K = 1 , NRTI
        WRITE(*,220) K, TI(K), (X(J,K),J=1,N)
      end do
  220 FORMAT(' ',I2,1X,F6.4,1P,5(2X,D12.5))
 3000 CONTINUE
c
c  Terminate.
c
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'musn_test():'
      write ( *, '(a)' ) '  Normal end of execution.'

      return
      end
      subroutine fdifn ( t, y, f )

c*********************************************************************72
c
cc fdifn evaluates the right hand side of the ODE.
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solution of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
c  Parameters:
c
c    Input, double precision T, the value of the independent variable.
c
c    Input, double precision Y(N), the dependent variable at T.
c
c    Output, double precision F(N), the right hand side at T, Y.
c
      implicit none

      double precision c1
      double precision c2
      double precision c3
      double precision c4
      double precision f(5)
      double precision t
      double precision y(5)

      c1 = 0.5D+00
      c2 = 0.9D+00
      c3 = 1000.0D+00
      c4 = 10.0D+00

      F(1) = c1 * Y(1) * ( Y(3) - Y(1) ) / Y(2)
      F(2) = - c1 * ( Y(3) - Y(1) )
      F(3) = ( c2 - c3 * ( Y(3) - Y(5) ) 
     &  - c1 * Y(3) * ( Y(3) - Y(1) ) ) / Y(4)
      F(4) = c1 * ( Y(3) - Y(1) )
      F(5) = c3 / c4 * ( Y(3) - Y(5) )

      return
      end
      subroutine x0t ( t, x )

c*********************************************************************72
c
cc X0T supplies the initial approximate solution to the nonlinear BVP.
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solutioni of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
c  Parameters:
c
c    Input, double precision T, the value of the independent variable.
c
c    Output, double precision X(N), the approximate solution at T.
c
      implicit none

      double precision t
      double precision X(5)

      X(1) = 1.0D+00
      X(2) = 1.0D+00
      X(4) = -10.0D+00
      X(3) = - 4.5D0*T*T + 8.91D+00 * T + 1.0D+00
      X(5) = - 4.5D0*T*T + 9.0D+00 * T + 0.91D+00

      return
      end
      subroutine G ( N, XA, XB, FG, DGA, DGB )

c*********************************************************************72
c
cc G evaluates the nonlinear boundary condition.
c
c  Modified:
c
c    09 January 2013
c
c  Author:
c
c    Robert Mattheij, G.W.M. Staarink.
c
c  Reference:
c
c    Uri Ascher, Robert Mattheij, Robert Russell,
c    Numerical Solutioni of Boundary Value Problems for 
c    Ordinary Differential Equations,
c    Prentice Hall, 1988,
c    ISBN: 0-13-627266-5,
c    LC: QA379.A83.
c
c  Parameters:
c
c    Input, integer N, the number of equations.
c
c    Input, double precision XA(N), XB(N), the value of the dependent
c    variables at A and B.
c
c    Output, double precision FG(N), the value of the boundary conditions.
c
c    Output, double precision DGA(N,N), DGB(N,N), the Jacobians of FG
c    with respect to XA(N) and XB(N) respectively.
c
      implicit none

      integer n

      double precision dga(n,n)
      double precision dgb(n,n)
      double precision fg(n)
      integer i
      integer j
      double precision xa(n)
      double precision xb(n)

      FG(1) = XA(1) - 1.0D0
      FG(2) = XA(2) - 1.0D0
      FG(3) = XA(3) - 1.0D0
      FG(4) = XA(4) + 10.0D0
      FG(5) = XB(3) - XB(5)

      DO I = 1, N
        DO J = 1, N
          DGA(I,J) = 0.0D0
          DGB(I,J) = 0.0D0
        end do
      end do

      DGA(1,1) = 1.0D0
      DGA(2,2) = 1.0D0
      DGA(3,3) = 1.0D0
      DGA(4,4) = 1.0D0
      DGB(5,3) = 1.0D0
      DGB(5,5) = -1.0D0

      RETURN
      END
