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  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
c   Make up output a file name that means something
c   Combining 'odeout-' with the step size
write(xstep,'(f5.3)') h
open(10,file=filename)
write(10,2010)
c
c     Get A number of integration steps
c
write(*,'(a)',advance='no') 'Number of integration steps: '
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.
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
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
yab=ynew
c
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
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.
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