subroutine get_unit ( iunit )

c*********************************************************************72
c
cc get_unit() returns a free FORTRAN unit number.
c
c  Discussion:
c
c    A "free" FORTRAN unit number is an integer between 1 and 99 which
c    is not currently associated with an I/O device.  A free FORTRAN unit
c    number is needed in order to open a file with the OPEN command.
c
c    If IUNIT = 0, then no free FORTRAN unit could be found, although
c    all 99 units were checked (except for units 5, 6 and 9, which
c    are commonly reserved for console I/O).
c
c    Otherwise, IUNIT is an integer between 1 and 99, representing a
c    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
c    are special, and will never return those values.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    15 January 2008
c
c  Author:
c
c    John Burkardt
c
c  Output:
c
c    integer IUNIT, the free unit number.
c
      implicit none

      integer i
      integer ios
      integer iunit
      logical lopen

      iunit = 0

      do i = 1, 99

        if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

          inquire ( unit = i, opened = lopen, iostat = ios )

          if ( ios == 0 ) then
            if ( .not. lopen ) then
              iunit = i
              return
            end if
          end if

        end if

      end do

      return
      end
      function humps_antideriv ( x )

c*********************************************************************72
c
cc humps_antideriv() evaluates the antiderivative of the humps function.
c
c  Discussion:
c
c    y = 1.0 / ( ( x - 0.3 )^2 + 0.01 )
c      + 1.0 / ( ( x - 0.9 )^2 + 0.04 )
c      - 6.0
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 August 2019
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision x: the argument.
c
c  Output:
c
c    double precision humps_antideriv: the value of the antiderivative at x.
c
      implicit none

      double precision humps_antideriv
      double precision x
      double precision ya

      ya = ( 1.0D+00 / 0.1D+00 ) * atan ( ( x - 0.3D+00 ) / 0.1D+00 )   
     &   + ( 1.0D+00 / 0.2D+00 ) * atan ( ( x - 0.9D+00 ) / 0.2D+00 )   
     &   - 6.0D+00 * x

      humps_antideriv = ya

      return
      end
      function humps_deriv ( x )

c*********************************************************************72
c
cc humps_deriv() evaluates the derivative of the humps function.
c
c  Discussion:
c
c    y = 1.0 / ( ( x - 0.3 )^2 + 0.01 )
c      + 1.0 / ( ( x - 0.9 )^2 + 0.04 )
c      - 6.0
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 August 2019
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision x: the argument.
c
c  Output:
c
c    double precision humps_deriv: the value of the derivative at x.
c
      implicit none

      double precision humps_deriv
      double precision x
      double precision yp

      yp = - 2.0D+00 * ( x - 0.3D+00 ) 
     &    / ( ( x - 0.3D+00 )**2 + 0.01D+00 )**2        
     &     - 2.0D+00 * ( x - 0.9D+00 ) 
     &    / ( ( x - 0.9D+00 )**2 + 0.04D+00 )**2

      humps_deriv = yp

      return
      end
      function humps_deriv2 ( x )

c*********************************************************************72
c
cc humps_deriv2() evaluates the second derivative of the humps function.
c
c  Discussion:
c
c    y = 1.0 / ( ( x - 0.3 )^2 + 0.01 )
c      + 1.0 / ( ( x - 0.9 )^2 + 0.04 )
c      - 6.0
c
c    ypp = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )^2 + 0.01 )^2
c          - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )^2 + 0.04 )^2
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 August 2019
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision x: the argument.
c
c  Output:
c
c    double precision humps_deriv2: the value of the second derivative at x.
c
      implicit none

      double precision humps_deriv2
      double precision u1
      double precision u1p
      double precision u2
      double precision u2p
      double precision v1
      double precision v1p
      double precision v2
      double precision v2p
      double precision x
      double precision ypp

      u1 = - 2.0D+00 * ( x - 0.3D+00 )
      v1 = ( ( x - 0.3D+00 )**2 + 0.01D+00 )**2
      u2 = - 2.0D+00 * ( x - 0.9D+00 )
      v2 = ( ( x - 0.9D+00 )**2 + 0.04D+00 )**2

      u1p = - 2.0D+00
      v1p = 2.0D+00 * ( ( x - 0.3D+00 )**2 + 0.01D+00 ) 
     &  * 2.0D+00 * ( x - 0.3D+00 )
      u2p = - 2.0D+00
      v2p = 2.0D+00 * ( ( x - 0.9D+00 )**2 + 0.04D+00 ) 
     &  * 2.0D+00 * ( x - 0.9D+00 )

      ypp = ( u1p * v1 - u1 * v1p ) / v1 / v1
     &    + ( u2p * v2 - u2 * v2p ) / v2 / v2

      humps_deriv2 = ypp

      return
      end
      function humps_fun ( x )

c*********************************************************************72
c
cc humps_fun() evaluates the humps function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 August 2019
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision x: the evaluation point.
c
c  Output:
c
c    double precision humps_fun: the function value.
c
      implicit none

      double precision humps_fun
      double precision x
      double precision y

      y = 1.0D+00 / ( ( x - 0.3D+00 )**2 + 0.01D+00 )
     &  + 1.0D+00 / ( ( x - 0.9D+00 )**2 + 0.04D+00 )
     &  - 6.0D+00

      humps_fun = y

      return
      end
      function humps_ode ( x, y )

c*********************************************************************72
c
cc humps_ode() evaluates the derivative of the humps function for an ODE solver.
c
c  Discussion:
c
c    This verion of "humps_deriv" appends the input argument "y", as expected
c    by most ODE solving software.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 August 2019
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision x: the argument.
c
c    double precision y: the value of the dependent variable.
c
c  Output:
c
c    double precision humps_ode: the value of the derivative of the humps function.
c
      implicit none

      double precision humps_ode
      double precision x
      double precision y
      double precision yp

      yp = - 1.0D+00 / ( ( x - 0.3D+00 )**2 + 0.01D+00 )**2
     &     * 2.0D+00 * ( x - 0.3D+00 )
     &     - 1.0D+00 / ( ( x - 0.9D+00 )**2 + 0.04D+00 )**2
     &     * 2.0D+00 * ( x - 0.9D+00 )

      humps_ode = yp

      return
      end
      subroutine plot_xy ( n, x, y, prefix )

c*********************************************************************72
c
cc plot_xy() plots xy data.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    06 August 2019
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer n : the number of data points.
c
c    double precision x(n), y(n): the data points.
c
c    character ( len = * ) prefix: the prefix for the plot names.
c
      implicit none

      integer n

      character ( len = 255 ) command_filename
      integer command_unit
      character ( len = 255 ) data_filename
      integer data_unit
      integer i
      character ( len = 255 ) output_filename
      character ( len = * ) prefix
      character ( len = 255 ) prefix2
      double precision x(n)
      double precision y(n)
c
c  Create the data file.
c
      data_filename = trim ( prefix ) // '_data.txt'
      call get_unit ( data_unit )
      open ( unit = data_unit, file = data_filename, 
     &  status = 'replace' )
      do i = 1, n
        write ( data_unit, '(2x,g14.6,2x,g14.6)' ) x(i), y(i)
      end do
      close ( unit = data_unit )
      write ( *, '(a)' ) ''
      write ( *, '(a)' )     '  Created graphics data file "' 
     &  // trim ( data_filename ) // '".'
c
c  Plot the selected data.
c
      command_filename = trim ( prefix ) // '_commands.txt'
      call get_unit ( command_unit )
      open ( unit = command_unit, file = command_filename, 
     &  status = 'replace' )

      write ( command_unit, '(a)' ) '# ' // trim ( command_filename )
      write ( command_unit, '(a)' ) '#'
      write ( command_unit, '(a)' ) '# Usage:'
      write ( command_unit, '(a)' ) '#  gnuplot < ' 
     &  // trim ( command_filename )
      write ( command_unit, '(a)' ) '#'
      write ( command_unit, '(a)' ) 'set term png'
      write ( command_unit, '(a)' ) 'set nokey'
      output_filename = trim ( prefix ) // '.png'
      write ( command_unit, '(a)' )     'set output "' 
     &  // trim ( output_filename ) // '"'
      write ( command_unit, '(a)' ) 'set xlabel "<---X--->"'
      write ( command_unit, '(a)' ) 'set ylabel "<---Y(X)--->"'
      call s_escape_tex2 ( prefix, prefix2 )
      write ( command_unit, '(a)' ) 'set title "' 
     &  // trim ( prefix2 ) // '"'
      write ( command_unit, '(a)' ) 'set grid'
      write ( command_unit, '(a)' ) 'set style data lines'
      write ( command_unit, '(a)' ) 'plot "' // trim ( data_filename )  
     &   // '" using 1:2 lw 3 linecolor rgb "blue"'

      close ( unit = command_unit )
      write ( *, '(a)' )     '  Created graphics command file "' 
     &  // trim ( command_filename ) // '".'

      return
      end
      subroutine s_escape_tex2 ( s1, s2 )

c*********************************************************************72
c
cc s_escape_tex2() de-escapes TeX escape sequences.
c
c  Discussion:
c
c    In particular, every occurrence of the characters '\', '_',
c    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
c    '\{' and '\}'.  A TeX interpreter, on seeing these character
c    strings, is then likely to return the original characters.
c
c    In some cases, it seems that TWO backslashes are needed.
c    This version of the function provides them.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, character ( len = * ) S1, the string to be de-escaped.
c
c    Output, character ( len = * ) S2, a copy of the string,
c    modified to avoid TeX escapes.
c
      implicit none

      character ch
      character ( len = * ) s1
      integer s1_length
      integer s1_pos
      character ( len = * ) s2
      integer s2_pos

      s1_length = len_trim ( s1 )

      s1_pos = 0
      s2_pos = 0
      s2 = ' '

      do while ( s1_pos < s1_length )

        s1_pos = s1_pos + 1
        ch = s1(s1_pos:s1_pos)

        if ( ch == '\\' .or.
     &       ch == '_' .or.
     &       ch == '^' .or.
     &       ch == '{' .or.
     &       ch == '}' ) then

          s2_pos = s2_pos + 1
          s2(s2_pos:s2_pos) = '\\'
          s2_pos = s2_pos + 1
          s2(s2_pos:s2_pos) = '\\'

        end if

        s2_pos = s2_pos + 1
        s2(s2_pos:s2_pos) = ch

      end do

      return
      end