program ode
      implicit none
*
* we define our variables here
*      
* make x, y arrays so we can plot them using pgplot
*
      integer maxits
      parameter (maxits = 10000)
      real h, a, b, x, y, z
      real x0, y0, z0
      real xmin, xmax, ymin, ymax
      real xplot(maxits), yplot(maxits) 
      real xreal(maxits), yreal(maxits)
      integer i, n
*
* let's experiment with accuray by reading in different n
*
      n = 100
      write(6,*)'enter a value for n (MAX 1000): '  
      read(5,*)n

*
* setting our initial conditions 
* we are integrating from x = a to x = b
      
      a = 0.0
      b = 1.0

* set h, our increment
*
      h = (b - a)/float(n)
*
* Setting the initial conditions:  y(a) = 1.0, dy/dx(a) = z
*
      x0 = a
      y0 = 1.0
      z0 = 0.0

* for plotting

      xplot(1) = x0
      yplot(1) = y0
*
* subroutine to solve the second order ODE
* (you may choose to use Euler's method for comparison!)

      call rungekutta2(n, h, x0, y0, z0,   xplot, yplot)
c      call euler2(n, h, x0, y0, z0, xplot, yplot)
*
* now get 'real' solution
*
      call realfunc(a, b, x0, y0, xreal, yreal)

* we'll write out the data from the rungekutta subroutine here
*
      open(10,file='rk.dat',status='unknown')
      write(10,*)'#      x           y'
      do i = 1, n
         write(10,99)xplot(i),yplot(i)
      enddo
 99   format(1x,f10.5,1x,f10.5)  ! notice the format statement, 10characters with 5 numbers 
                                 ! after the decimal point
      close(10)
*
* and the data from the exact solution here
*
      open(12,file='exact.dat',status='unknown')
      write(12,*)'#      x           y'
      do i = 1, 1000
         write(12,99)xreal(i),yreal(i)  ! we can use same format statement
      enddo
      close(12)

      stop
      end

* using euler's method
*
      subroutine euler2(n, h, x0, y0, z0, xplot, yplot)
      implicit none
      real h, x0, y0, z0, xplot(1000), yplot(1000)
      real xl, yl, zl, x, y, z, f
      integer i, n
      x = x0
      y = y0
      z = z0
      do i = 1, n
         xl = x
         yl = y
         zl = z
         x = xl + h
         y = yl + h*zl
         z = zl + h*f(xl, yl, zl)
         xplot(i) = x
         yplot(i) = y
      enddo
      return
      end

* subroutine that solves the 2nd order ODEs using the Runge-Kutta method
*
      subroutine rungekutta2(n, h, x0, y0, z0,   xplot, yplot)
      implicit none
      real h, x0, y0, z0, xplot(1000), yplot(1000)
      real x, y, z, xl, yl, zl, k1, k2, k3, k4, f
      real m1, m2, m3, m4, g
      integer i, n
      x = x0
      y = y0
      z = z0
      do i = 1, n
         xl = x
         yl = y
         zl = z
         k1 = h*f(xl,yl,zl)
         m1 = h*g(xl,yl,zl)
         k2 = h*f((xl+0.5*h),(yl+k1*0.5),(zl+m1*0.5))
         m2 = h*g((xl+0.5*h),(yl+k1*0.5),(zl+m1*0.5))
         k3 = h*f((xl+0.5*h),(yl+k2*0.5),(zl+m2*0.5))
         m3 = h*g((xl+0.5*h),(yl+k2*0.5),(zl+m2*0.5))
         k4 = h*f((xl+h),(yl+k3),(zl+m3))
         m4 = h*g((xl+h),(yl+k3),(zl+m3))
         x = xl + h
         y = yl + 0.16666666*(k1 + 2.*k2 + 2.*k3 + k4)
         z = zl + 0.16666666*(m1 + 2.*m2 + 2.*m3 + m4)
c         write(10,*)x,y,z
         xplot(i) = x
         yplot(i) = y
      enddo
      return
      end

* subroutine which gives x,y arrays for analytical function
*
      subroutine realfunc(a, b, x0, y0,     x, y)
      implicit none
      real x0, y0, x(1000), y(1000), dx, a, b
      integer i
      x(1) = x0
      y(1) = y0
      dx = (b-a)/float(1000)
      do i = 2, 1000
         x(i) = i*dx
         y(i) = 0.25*(exp(2.0*x(i))*(2.0*x(i) - 1.0)) + 1.25
      enddo
      return
      end

*
* function for g = exp(2x) - 2z
*
      real function g(x,y,z)
      implicit none
      real x,y,z
      g = exp(2.0*x)+2.0*z
      return
      end
*
* function for z = dy/dx
*      
      real function f(x,y,z)
      implicit none
      real x,y,z
      f = z
      return
      end