FORTRANの場合

      PROGRAM shot
      DIMENSION r(2),v(2),a(2)
      CALL init(t,r,v,a,dt,ncalc,g,c)
      CALL prtout(t,r,v,a,0)
      CALL euler(t,r,v,a,0.0,1,g,c)
      CALL prtout(t,r,v,a,1)
      DO 10 i = 1,100
        CALL euler(t,r,v,a,dt,ncalc,g,c)
        CALL prtout(t,r,v,a,1)
        IF(r(2) .LE. 0.0) GOTO 999
   10 CONTINUE
  999 END
C
      SUBROUTINE init(t,r,v,a,dt,ncalc,g,c)
      DIMENSION r(2),v(2),a(2)
      g = 9.81
      t = 0
      prper = 0.2
      dt = 0.001
*      WRITE(*,*) 'プリントの間隔(s), 時間の刻み巾(s) = '
*      READ(*,*) prper,dt
      ncalc = prper/dt + 0.5
      height = 10
      v0 = 20
      th = 40
*      WRITE(*,*) '初期条件： 高さ(m), 速さ(m/s), 角度(degree) = '
*      READ(*,*) height,v0,th
      th = th/180*3.141592
      WRITE(*,*) '空気抵抗の係数 (1/m) = '
      READ(*,*) c
      v(1) = v0*cos(th)
      v(2) = v0*sin(th)
      r(1) = 0
      r(2) = height
      END
C
      SUBROUTINE prtout(t,r,v,a,n)
      DIMENSION r(2),v(2),a(2)
      IF(n.EQ.0) THEN
        WRITE(*,'(/2(1X,2A/))') 
     +   '   t       x       y      v       a   ',
     +   '    vx      vy      ax      ay  ',
     +   '  (s)     (m)     (m)   (m/s)  (m/s^2)',
     +   '  (m/s)   (m/s)  (m/s^2) (m/s^2)'
      ELSE
        va = sqrt(v(1)**2+v(2)**2)
        aa = sqrt(a(1)**2+a(2)**2)
        WRITE(*,'(1X,F5.2,8F8.2)') t,r,va,aa,v,a
      ENDIF
      END
C
      SUBROUTINE euler(t,r,v,a,dt,ncalc,g,c)
      DIMENSION r(2),v(2),a(2)
      DO 10 i = 1, ncalc
        r(1) = r(1) + v(1)*dt
        r(2) = r(2) + v(2)*dt
        vv = sqrt(v(1)**2 + v(2)**2)
        a(1) = - c*vv*v(1)
        a(2) = - g -c*vv*v(2)
        v(1) = v(1) + a(1)*dt
        v(2) = v(2) + a(2)*dt
        t = t + dt
        IF(r(2) .LE. 0.0) RETURN
   10 CONTINUE
      END