c codeint.f c c
c
      module precision
c
c
         integer, parameter :: rp=selected_real_kind(13,300)
c
      end module precision
      program odeint
      use precision
      implicit none
c
c  Demonstration of Runge-Kutta, Adams-Bashford, and Adams-Moulton for
c  solution of an Ordinary Differential Equation in the form:
c
c     d y
c     ---  =  f ( x , y )
c     d x
c
c  The arrays below allow storage of the results for the four most recent
c  time steps.  For example y(0) contains the value of y at the beginning
c  of the current time step, y(1) contains the value of y at the beginning
c  or the previous time step, y(2) contains the value of y 2 steps back, etc.
c
      real (rp) f(0:3),y(0:3),fa(0:3),ya(0:3),fimp(0:3),yt(0:3),ft(0:3)
      real (rp) h, ynew, yrk, yab, yanew, xnew,  dfdy, xn, xk1,
     &          xk2, xk3, xk4, yam, yimpn, yimp, gy, dgdy, dely, yabt,
     &          fnew, func, ftnew
      integer nsteps, istep, i, it
      character*12 filename
      character*5 xstep
c   Read in the integration step size
c 
    1 write(*,'(a)',advance='no') 'Provide Integration Step Size: '
c  
      read (*,*)h
c   Make up output a file name that means something
c   Combining 'odeout-' with the step size
      write(xstep,'(f5.3)') h
      filename = 'odeout-'//adjustl(xstep)
      open(10,file=filename)
      write(10,2010)
 2010 format(4x,'x',10x,'correct',5x,'Runge',7x,'Adams',7x,'Adams',4x,
     $  'Implicit',/15x,'answer',6x,'Kutta',7x,'Bashford',4x,'Moulton',
     $  2x,'Adams-Moultn',/)
c
c     Get A number of integration steps
c
      write(*,'(a)',advance='no') 'Number of integration steps: '
      read(*,*) nsteps
c
c    Take at least 5 steps
c
      nsteps=max(nsteps,5)
c
c    Here We hardwire the initial conditions
c
      ynew=0.
      xnew=0.
      call answer(xnew,yanew)
c
c   Initialize values of previous timesteps
c
      do 10 i=1,3
      y(i)=0.
      ya(i)=0.
      yt(i)=0.
      f(i)=0.
      fimp(i)=0.
      ft(i)=0.
   10 fa(i)=0.
c
c   Do Three steps of Runge Kutta
c
      do 40 istep=1,3
c
c   Shuffle past data
c
      	do 30 i=3,1,-1
           yt(i)=yt(i-1)
           ya(i)=ya(i-1)
      	   y(i)=y(i-1)
           f(i)=f(i-1)
           ft(i)=ft(i-1)
           fimp(i)=fimp(i-1)
   30      continue
        y(0)=ynew
        yt(0)=ynew
        ya(0)=yanew
      	xn= xnew
c
c    Make evaluations for this time step
c
        f(0)=func(xn,y(0),dfdy)
        ft(0)=f(0)
c
c   Initialize numbers needed for the Implicit Adams-Moulton method.
c   Note that I am cheating here  fimp should be
c   loaded with the results of an implicit Runge-Kutta
c   or explicit formulation for variable step size using
c   smaller steps for initialization
c
        fimp(0)=func(xn,ya(0),dfdy)
c
c   Actual  Runge-Kutta step
c
      	xk1=h*func(xn,y(0),dfdy)
      	xk2=h*func(xn+.5*h,y(0)+.5*xk1,dfdy)
      	xk3=h*func(xn+.5*h,y(0)+.5*xk2,dfdy)
      	xk4=h*func(xn+h,y(0)+xk3,dfdy)
c
c   Values at the new time level (end of time step)
c   If you look at this carefully it is a close relative of a
c   Simpson's Rule integration
c
        ynew=y(0)+(xk1+2.*xk2+2.*xk3+xk4)/6.
        xnew=xn+h
c
c   Get the actual answer for comparison
c
        call answer(xnew,yanew)
        write(10,2000)xnew,yanew,ynew
   40 continue
 2000 format(1p,6e12.5)
c
c    Now that we have starter values for Adams-Bashford and
c    Adams-Moulton, Finish integration with all methods
c
c    Current y for Runge-Kutta
c
      yrk=ynew
c
c    Current y for Adams-Bashford
c
      yab=ynew
c
c    Current y for Adams-Moulton
c
      yam=ynew
c
c    Current y for Implicit Adams-Moulton
c
      yimpn=yanew
c
c
c
      do 80 istep=4,nsteps
c
c   Set up for next time step by shifting stored values of previous time
c   steps to keep only the results of the three previous steps
c
      	do 60 i=3,1,-1
           ya(i)=ya(i-1)
           yt(i)=yt(i-1)
      	   y(i)=y(i-1)
           f(i)=f(i-1)
           ft(i)=ft(i-1)
           fimp(i)=fimp(i-1)
   60      continue
      	xn= xnew
        y(0)=yam
        yt(0)=yab
        f(0)=func(xn,y(0),dfdy)
        ft(0)=func(xn,yt(0),dfdy)
c
c   The Implicit Adams-Moulton requires a Newton Iteration (DO loop 70)
c   to solve a non-linear equation in the new-time value of y
c
        yimp=yimpn
        fimp(0)=func(xn,yimp,dfdy)
        xnew=xn+h
        do 70 it=1,10
          fnew=func(xnew,yimpn,dfdy)
          gy=yimp-yimpn+(9.*fnew+19.*fimp(0)-5.*fimp(1)+fimp(2))*h/24.
          dgdy=9.*h/24.*dfdy-1.
          dely=-gy/dgdy
          yimpn=yimpn+dely
c         
	  if(abs(dely/max(yimp,yimpn)).lt.1.e-9) go to 72
c           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
   70   continue
        write(6,*) 'Implicit iteration failed'
        stop
c
c   Standard Runge-Kutta continues
c
   72 	xk1=h*func(xn,yrk,dfdy)
      	xk2=h*func(xn+.5*h,yrk +.5*xk1,dfdy)
      	xk3=h*func(xn+.5*h,yrk +.5*xk2,dfdy)
      	xk4=h*func(xn+h,yrk +xk3,dfdy)
        yrk =yrk +(xk1+2.*xk2+2.*xk3+xk4)/6.
c
c   Do a straight Adams-Bashford integration
c
        yab=yam+(55.*ft(0)-59.*ft(1)+37.*ft(2)-9.*ft(3))*h/24.
c
c   Do an Adams-Bashford predictor for Adams-Moulton
c
        yabt=yam+(55.*f(0)-59.*f(1)+37.*f(2)-9.*f(3))*h/24.
        ftnew=func(xnew,yabt,dfdy)
c
c   Now apply and Adams-Moulton correction step
c
        yam=yam+(9.*ftnew+19.*f(0)-5.*f(1)+f(2))*h/24.
        call answer(xnew,yanew)
        write(10,2000)xnew,yanew,yrk,yab,yam,yimpn
   80 continue
c
      close(10)
      stop
      end
c
      function func(x,y,dfdy)
      use precision
      real (rp) x,y,dfdy,func
c
c   A specific function f for the right hand side of the ODE
c
         func=1.-y**2
c
c   Limit the value to keep the code running when one numerical solution
c   method goes nuts
c
         func=min(1.e10_rp,max(func,-1.e10_rp))
         dfdy=-2*y
      return
c
      end function func
c
      subroutine answer(x,y)
      use precision
      real (rp) x,y,e2x
c
c  For the function in "func" this is the analytic solution to the ODE
c
      e2x=exp(2*x)
c
      y=(e2x-1.)/(1.+e2x)
      return
      end
c
c c