c     This program solves the initial value problem:
c
c     dx/dt=f(x,t)
c
c     using the Runga Kutta order 4 method
c
c
c     In this case, the problem solved is:
c
c     dx/dt=1+x
c     x(1)=1
c
c     where in this case the function is time independent.
c     The numerical solution is compared to the true solution
c
c     xtrue=2exp(t)-1
c
c     Variables used are:
c
c     x ---  x computed by Runka Kutta
c     xtrue --- exact solution
c     t ---  time variable
c     f(x) ---  function (equal to 1+x)
c     dt --- time step of 0.1
c***************************************************

      parameter(ni=4)
      real k1,k2,k3,k4
c***************************************************
c Define the function at beginning of program.
c FORTRAN understands this is a function statement
c if placed at beginning.
c***************************************************
      f(x,t)=2.+(x-t-1.)**2
c***************************************************
c Set up time step and initial condition
c***************************************************
      delt=0.1
      x=2
      t=1
c***************************************************
c Output initial values. This will nee to be changed if function changed.
c***************************************************
      xtrue=1.+t+tan(t-1)
      write(*,900)
  900 format('t',5x,'x',6x,'xtrue',2x,'xtrue-x')
      write(*,1000) t,x,xtrue,(xtrue-x)
 1000 format(f3.1,2x,f7.5,1x,f7.5,1x,f7.5)
c***************************************************
c Perform RK4
c***************************************************
      do i=1,ni
         k1=delt*f(x,t)
         print*,k1
         k2=delt*f(x+k1/2,t+delt/2)
         print*,k2
         k3=delt*f(x+k2/2,t+delt/2)
         print*,k3
         k4=delt*f(x+k3,t+delt)
         print*,k4
         x=x+(k1+2.*k2+2.*k3+k4)/6.
         t=t+delt
c***************************************************
c Compare to analytical solution. Change or remove if different function used.
c***************************************************
         xtrue=1.+t+tan(t-1)
         write(*,1000) t,x,xtrue,(xtrue-x)
      enddo
      end