!
MODULE Harmonic
      USE ODE, ONLY: wp
      IMPLICIT NONE
      PUBLIC
      INTEGER, PARAMETER :: n_eq = 2
      REAL :: k = 1.0
CONTAINS
      SUBROUTINE kx (x, y, dy)
         INTEGER, PARAMETER :: wp = KIND (0.0D0)
         REAL (KIND=wp), INTENT (IN) :: x
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (IN) :: y
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (OUT) :: dy
         dy = (/ - k * y (2), y (1) /)
      END SUBROUTINE kx
END MODULE Harmonic
!
PROGRAM ODE_test
      USE Harmonic
      USE ODE
      USE SimpleGraphics
!
      INTEGER, PARAMETER :: sp = KIND (0.0E0), n_pts = 50
      REAL (KIND=wp), DIMENSION (n_eq, n_pts) :: x
      REAL (KIND=wp), DIMENSION (n_pts) :: t
      REAL (KIND=wp) :: t_min = 0.0_wp, t_max = 10.0_wp
      INTEGER :: i
      CHARACTER (50), DIMENSION (2) :: title ! The title of the plot
!
      CALL ODESolve (f=kx, n_eq=n_eq, x_range= (/ t_min, t_max /), n_points=n_pts, y0= (/ 0.0_wp, 1.0_wp /), x=t, y=x, relerr=1D-5)
!
  ! write(*,*) t
  ! write(*,*) "Velocity: ", x(1,:)
  ! write(*,*) "Position: ", x(2,:)
!
  ! This is an example of using Plot2D plot:
      title (1) = "x(t) for an anharmonic oscillator" ! First line of the title
      title (2) = "F(x)=-k*x+k/2*x^3" ! Second line of title
      CALL InitGraphics (file="Harmonic.ps", file_type="CONS", plot_title=title, x_label="t", y_label="x(t)")
      CALL Plot2D (x=REAL(t, sp), y=REAL(x(1, :), sp), plot_spec="L-R", new_plot=.TRUE.)! x=f(t)
      CALL EndGraphics ()
!
END PROGRAM ODE_test
!