      program prog083
      .....
      do i = 1, 89
        theta = dble(i)
        call solvexf(v0,theta,ga,dt,xf,tf)
        write(*,100) theta, xf, tf
      end do
      ...
      end
c **********************************************************************
      subroutine solvexf(v0,theta,ga,dt,xf,tf)
c **********************************************************************
      .....
      pi = .....
       t = 0.0d0
       x = .....
      .....
 10   continue
        call onestep(x,y,vx,vy,t,dt,ga)
        .....
      if( y .ge. 0.0d0 ) goto 10
      call solvey0(x,y,vx,vy,t,ga)
      xf = x
      tf = t
c
      end
