*************************************************************************

      subroutine lobbie(x,y,z,ts,tf,step,etol,nxy,nz,ns,ni,nst,ncf,fun)

*************************************************************************
*                                                                       *
*     The integrator Lobbie III numerically solves mixed systems of     *
*     ordinary differential equations of second and first orders        *
*     x''=f(t,x,x',z) & z'=g(t,x,x',z). The integrator is a software    *
*     implimentation of a collocation Lobatto method based on           *
*     Newton's polynomial interpolants for the functions f & g.         *
*     The theoretical background of the method is presented in          *
*     Avdyushev V.A. New Collocation Integrator for Solving Dynamic     *
*     Problems // Russ. Phys. J. Vol. 63. P. 1977–1988 (2021).          *
*                                                                       *
*     Input parameters :                                                *
*                                                                       *
*     x,y,z  initial values for x,x',z at ts                            *
*     ts     initial value for t                                        *
*     tf     final value for t                                          *
*     step   start stepsize                                             *
*     etol   local tolerances                                           *
*     nxy    dimension of x & y                                         *
*     nz     dimension of z                                             *
*     ns     number of nodes (2 <= ns <= 17)                            *
*     ni     number of iterations                                       *
*     fun    name of subroutine for f & g                               *
*                                                                       *
*     Notes :                                                           *
*                                                                       *
*     1. The order of the integration scheme is p = 2 ns - 2.           *
*     2. If step = 0, the start stepsize is selected by                 *
*        the integrator as for ns = 2 (p = 2).                          *
*     3. If all etol = 0, integration runs with the constant            *
*        start stepsize.                                                *
*     4. If ni < 1, iterative process runs until convergence by x       *
*        (up to 25 iterations).                                         *
*                                                                       *
*     Output parameters :                                               *
*                                                                       *
*     x,y,z  final values for x,x',z at tf                              *
*     step   penultimate stepsize                                       *
*     nst    number of steps                                            *
*     ncf    number of calls of fun                                     *
*                                                                       *
*     Subroutine for f & g :                                            *
*                                                                       *
*     fun(t,x,y,z,fg)                                                   *
*                                                                       *
*     Input parameters :                                                *
*                                                                       *
*     x,y,z  values for x,x',z at t                                     *
*                                                                       *
*     Output parameters :                                               *
*                                                                       *
*     fg     values for f & g (dim fg = nxy + nz)                       *
*                                                                       *
*************************************************************************

      implicit real*8(a-h,o-z); logical ls,vs
      parameter (ni0=25,sgm=3.16d0,eps=1d-14)
      dimension a(ns,nxy+nz),b(ns,ns),c(ns),d(ns,ns),e(ns,ns)
      dimension u(nxy),v(nxy),w(nz),x(nxy),y(nxy),z(nz)
      dimension ax(ns,nxy+nz),da(ns,nxy+nz),us(nxy),g(ns+1)
      dimension fg(nxy+nz),fx(nxy),etol(nxy+nz),ietol(nxy+nz)
      common/spacing/spacing(56)

      external fun

      if(ts==tf) return

      if(ns<2.or.ns>17) stop 'Inadmissible Order!'
      
      ! Spacing

      nr=ns/2-1; l=((ns-1)/2-1)*(ns/2-1)
      c(1)=0d0; c(ns)=1d0; if(ns>2) c(nr+2)=0.5d0
      c(2:nr+1)=spacing(l+1:l+nr)
      c(ns-1:ns-nr:-1)=1d0-spacing(l+1:l+nr)

      ! Node Differences
      
      d=0d0
      do i=1,ns
      do j=1,ns
      d(i,j)=c(i)-c(j)
      end do
      end do

      ! Constants

      do i=1,ns

      prd=1d0
      do k=1,ns+1
      prd=prd*c(i)/dfloat(k) 
      g(k)=prd
      end do
      b(i,1)=g(1)
      e(i,1)=g(2)

      do j=2,ns
      do k=1,ns-j+2
      g(k)=d(i,j-1)*g(k)-dfloat(k)*g(k+1)
      end do
      b(i,j)=g(1)
      e(i,j)=g(2)
      end do

      end do

      ! Initialization
      
      a=0d0; ax=0d0; da=0d0; us=0d0; nzl=nxy+1; nzr=nxy+nz
      t=ts; h=dabs(step)*(tf-ts)/dabs(tf-ts); r=1d0
      ncf=1; ncf0=ni0*(ns-1); nst=0; if(ni<=0) ni=ni0

      ! Variable Step-Size Mode

      netol=0
      do i=1,nzr
      if(etol(i)/=0d0) then
      netol=netol+1; ietol(netol)=i
      end if
      end do
      
      vs=netol/=0; deg=1d0/dfloat(ns)
      rmx=sgm**deg; rmn=1d0/rmx

      call fun(t,x,y,z,a(1,:))

      ! Start Step-Size

      if(h==0d0) then
      if(.not.vs) stop 'Zero Step!'
      eta=dsqrt(eps); ecal=0d0
      do while(ecal==0d0)
      eta=eta*10d0; u=x+eta*y+eta**2*a(1,1:nxy)/2d0
      v=y+eta*a(1,1:nxy); w=z+eta*a(1,nzl:nzr)
      call fun(t+eta,u,v,w,fg); ncf=ncf+1
      do i=1,netol
      j=ietol(i)
      ecal=ecal+((fg(j)-a(1,j))/etol(j))**2
      end do
      end do
      h=dsqrt(2d0*eta/dsqrt(ecal))
      h=h*(tf-ts)/dabs(tf-ts)
      end if

    1 h=r*h; h2=h**2; ls=.false.

      ! First & Last Step

      if((tf-t)/h<1d0) then
      h=tf-t; h2=h**2; ls=.true.
      end if

      ! Integration
      
      do

      ! Iterations

    2 do it=1,ni

      do i=2,ns
      
      ! Collocation Solution

      fx=0d0
      do j=1,ns
      fx=fx+e(i,j)*a(j,1:nxy)
      end do
      u=x+h*c(i)*y+h2*fx

      fg=0d0
      do j=1,ns
      fg=fg+b(i,j)*a(j,:)
      end do
      v=y+h*fg(  1:nxy)
      w=z+h*fg(nzl:nzr)

      ! Function

      call fun(t+h*c(i),u,v,w,fg); ncf=ncf+1

      ! Divided Difference

      do j=1,i-1
      fg=(fg-a(j,:))/d(i,j)
      end do
      a(i,:)=fg

      end do

      ! Convergence Check

      if(all(dabs(u-us)<dabs(eps*u))) exit   
      us=u

      end do

      ! Forced Convergence
      
      if(nst==0.and.ncf<ncf0.and.it>ni) goto 2

      ! Scale Factor

      if(vs) then

      ecal=0d0
      do i=1,netol
      j=ietol(i)
      ecal=ecal+(a(ns,j)/etol(j))**2
      end do
      if(ecal/=0d0) then
      ecal=dabs(h)*dsqrt(ecal)*deg
      r=(1d0/ecal)**deg
      else
      r=rmx
      end if

      ! First Step-Size Adjustment

      if(nst==0) then
      if(rmx<r.and.ls  ) goto 3
      if(rmx<r.or.r<rmn) goto 1
      end if

      if(r>rmx) r=rmx ! Damping

      end if
      
      ! Solution

    3 fx=0d0
      do j=1,ns
      fx=fx+e(ns,j)*a(j,1:nxy)
      end do
      x=x+h*y+h2*fx

      fg=0d0
      do j=1,ns
      fg=fg+b(ns,j)*a(j,:)
      end do
      y=y+h*fg(  1:nxy)
      z=z+h*fg(nzl:nzr)
      
      t=t+h; nst=nst+1

      if(ls) exit ! Exit

      h=r*h; h2=h**2; us=0d0

      ! Last Step

      if((tf-t)/h<1d0) then
      r=r*(tf-t)/h; step=dabs(h) 
      h=tf-t; h2=h**2; ls=.true.
      end if

      ! Extrapolation
      
      do i=1,ns
      
      fg=0d0
      do j=ns,1,-1
      fg=a(j,:)+(1d0+r*c(i)-c(j))*fg
      end do

      do j=1,i-1
      fg=(fg-ax(j,:))/d(i,j)
      end do
      
      ! Correction

      if(nst/=1) da(i,:)=a(i,:)-ax(i,:)
      
      ax(i,:)=fg
      
      end do

      a=ax+da ! Predictor

      end do

      end subroutine lobbie
************************************************************************
      block data
      real*8, dimension(56) :: spacing
      common/spacing/spacing
      data spacing/
     *0.27639320225002103035908263312687,
     *0.17267316464601142810085377187657,
     *0.11747233803526765357449851302033,
     *0.35738424175967745184292450297956,
     *0.08488805186071653506398389301626,
     *0.26557560326464289309811405904562,
     *0.06412992574519669233127711938966,
     *0.20414990928342884892774463430102,
     *0.39535039104876056561567136982732,
     *0.05012100229426992134382737779083,
     *0.16140686024463112327705728645432,
     *0.31844126808691092064462396564567,
     *0.04023304591677059308553366958883,
     *0.13061306744724746249844691257008,
     *0.26103752509477775216941245363437,
     *0.41736052116680648768689011702091,
     *0.03299928479597043283386293195030,
     *0.10775826316842779068879109194577,
     *0.21738233650189749676451801526112,
     *0.35212093220653030428404424222047,
     *0.02755036388855888829620993084839,
     *0.09036033917799666082567920914154,
     *0.18356192348406966116879757277817,
     *0.30023452951732553386782510421652,
     *0.43172353357253622256796907213015,
     *0.02334507667891804405154726762227,
     *0.07682621767406384156703719645062,
     *0.15690576545912128696362048021682,
     *0.25854508945433189912653138318153,
     *0.37535653494688000371566314981289,
     *0.02003247736636954932244991899228,
     *0.06609947308482637449988989854586,
     *0.13556570045433692970766379973955,
     *0.22468029853567647234168864707046,
     *0.32863799332864357747804829817916,
     *0.44183406555814806617061164513192,
     *0.01737703674808071360207430396519,
     *0.05745897788851185058729918425888,
     *0.11824015502409239964794076201185,
     *0.19687339726507714443823503068163,
     *0.28968097264316375953905153063071,
     *0.39232302231810288088716027686354,
     *0.01521597686489103352387863081627,
     *0.05039973345326395350268586924007,
     *0.10399585406909246803445586451842,
     *0.17380564855875345526605839017970,
     *0.25697028905643119410905460707656,
     *0.35008476554961839595082327263885,
     *0.44933686323902527607848349747704,
     *0.01343391168429084292151024906313,
     *0.04456000204221320218809874680113,
     *0.09215187438911484644662472338123,
     *0.15448550968615764730254032131377,
     *0.22930730033494923043813329624797,
     *0.31391278321726147904638265963237,
     *0.40524401324084130584786849262344/
      end block data
*************************************************************************

      subroutine fun(t,x,y,z,fg)

*************************************************************************
*                                                                       *
*     This subroutine calculates the functions f & g of                 *
*     the differential equations x''=f(x) & z'=g(z) of the two-body     *
*     problem (dim x = dim f = 2; dim z = dim g = 4; dim fg = 6).       *
*                                                                       *
*************************************************************************

      implicit real*8(a-h,o-z)
      dimension x(2),y(2),z(4),fg(6)
      rx=dsqrt(x(1)**2+x(2)**2)
      rz=dsqrt(z(1)**2+z(2)**2)
      fg(1:2)=-x(1:2)/rx**3
      fg(3:4)= z(3:4)
      fg(5:6)=-z(1:2)/rz**3
      end subroutine fun
