BASICの場合

'PROGRAM shot
DIM 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, 1, g, c)
CALL Prtout(t, r(), v(), a(), 1)
LET i = 0
DO
  LET i = i + 1
  CALL Euler(t, r(), v(), a(), dt, ncalc, g, c)
  CALL Prtout(t, r(), v(), a(), 1)
LOOP UNTIL i > 100 OR r(2) < 0
END

SUB Euler (t, r(), v(), a(), dt, ncalc, g, c)
FOR i = 1 TO ncalc
  LET r(1) = r(1) + v(1) * dt
  LET r(2) = r(2) + v(2) * dt
  LET vv = SQR(v(1) ^ 2 + v(2) ^ 2)
  LET a(1) = -c * vv * v(1)
  LET a(2) = -g - c * vv * v(2)
  LET v(1) = v(1) + a(1) * dt
  LET v(2) = v(2) + a(2) * dt
  LET t = t + dt
IF r(2) <= 0 THEN EXIT SUB
NEXT i
END SUB

SUB Init (t, r(), v(), a(), dt, ncalc, g, c)
LET g = 9.81
LET t = 0
LET prper = .2
LET dt = .001
'INPUT "プリントの間隔(s), 時間の刻み巾(s) = "; prper,dt
LET ncalc = prper / dt
LET height = 10
LET v0 = 20
LET th = 40
'INPUT "初期条件： 高さ(m), 速さ(m/s), 角度(degree) = "; height,v0,th
LET th = th / 180 * 3.141592
INPUT "空気抵抗の係数 (1/m) = "; c
LET v(1) = v0 * COS(th)
LET v(2) = v0 * SIN(th)
LET r(1) = 0
LET r(2) = height
END SUB

SUB Prtout (t, r(), v(), a(), n)
IF n = 0 THEN
  PRINT
  PRINT "   t       x       y      v       a   ";
  PRINT "    vx      vy      ax      ay  "
  PRINT "  (s)     (m)     (m)   (m/s)  (m/s^2)";
  PRINT "  (m/s)   (m/s)  (m/s^2) (m/s^2)"
  PRINT
ELSE
  LET va = SQR(v(1) ^ 2 + v(2) ^ 2)
  LET aa = SQR(a(1) ^ 2 + a(2) ^ 2)
  PRINT USING "##.##"; t;
  PRINT USING " ####.##"; r(1); r(2); va; aa; v(1); v(2); a(1); a(2)
END IF
END SUB