Written in QBasic
download now

DECLARE SUB getvvector (pos.x#, pos.y#, masse.x#, masse.y#, masse.grav#, t!, vv.x#, vv.y#)
SCREEN 12
pos.x# = 300
pos.y# = 100
masse.x#(1) = 200
masse.y#(1) = 200
masse.grav#(1) = 4000
masse.x#(2) = 400
masse.y#(2) = 150
masse.grav#(2) = 4000
masse.x#(3) = 300
masse.y#(3) = 300
masse.grav#(3) = 4000
masse.x#(4) = 100
masse.y#(4) = 300
masse.grav#(4) = 4000
CIRCLE (masse.x#(1), masse.y#(1)), 1, 14
CIRCLE (masse.x#(2), masse.y#(2)), 1, 14
CIRCLE (masse.x#(3), masse.y#(3)), 1, 14
CIRCLE (masse.x#(4), masse.y#(4)), 1, 14
vk.x# = 10: vk.y# = 0: drawptstep# = 100: t = .0005
FOR i = 0 TO 60000 STEP t
k% = 1
DO UNTIL masse.x#(k%) = 0
CALL getvvector(pos.x#, pos.y#, masse.x#(k%), masse.y#(k%), masse.grav#(k%), t, vv.x#, vv.y#)
vk.x# = vk.x# + vv.x#
vk.y# = vk.y# + vv.y#
k% = k% + 1
LOOP
pos.x# = pos.x# + .5 * vk.x# * t
pos.y# = pos.y# + .5 * vk.y# * t
IF ABS(i / t) MOD drawptstep# = 0 THEN
IF 0 < pos.x# AND pos.x# < 700 AND 0 < pos.y# AND pos.y# < 500 THEN PSET (pos.x#, pos.y#)
END IF
NEXT i

SUB getvvector (pos.x#, pos.y#, masse.x#, masse.y#, masse.grav#, t, vv.x#, vv.y#)
d# = SQR((pos.x# - masse.x#) ^ 2 + (pos.y# - masse.y#) ^ 2)
am.l# = masse.grav# / d# ^ 2
am.x# = -(am.l# / d#) * (pos.x# - masse.x#)
am.y# = -(am.l# / d#) * (pos.y# - masse.y#)
vv.x# = am.x# * t
vv.y# = am.y# * t
END SUB