At the risk of being downvoted and rebuffed as waffling, I believe it would be useful for the OP and anyone encountering this question as a result of a search on similar terms, to look at a quick summary on the simplest way to numerically simulate the effect of forces on a point-like object.
(Why point-like? So that we can ignore any internal structure.)
First thing to realize, is that in numerical simulations, time is discretized to short intervals. Often, the intervals are of fixed (game-time) duration, to simplify the math involved. The key thing, however, is to realize that each step does not correspond to one instance of simulated continuous time, but to an interval. In mathematical terms, we need to integrate the instantaneous effects over the interval, to calculate (approximate) the effects during that interval.
Let's assume that for each object, at the start of each time step, we know their location/position vector $\vec{x}$, their velocity vector $\vec{v}$, and the sum of all forces that apply to it for the duration of the time step (about to be simulated), $\vec{F}$. If you have more than one force, you simply sum their vector components together. Let's use $\tau$ to describe the duration of the time step.
The simulated object has momentum $\vec{p} = m \vec{v}$, where $m$ is the mass of the object. If the mass of the object and the force acting on the object stay constant during the time step, then according to Newton's second law of motion,
$$m \vec{a} = \vec{F} \iff \vec{a} = \frac{1}{m} \vec{F}$$
where $\vec{a}$ is the acceleration during the time step. Noting that
$$\vec{v}(t) = \int_{0}^{t} \vec{a}(\tau) d \tau = \vec{v}(0) + \frac{t}{m} \vec{F}$$
and that
$$\vec{x}(t) = \int_{0}^{t} \vec{v}(\tau) d \tau = \vec{x}(0) + \vec{v}(0) t + \vec{F} \frac{t^2}{2 m}$$
In discrete time units, this means that
$$\vec{x} \leftarrow \vec{x} + \vec{v} \tau + \vec{F} \frac{\tau^2}{2 m}$$
and
$$\vec{v} \leftarrow \vec{v} + \frac{\tau}{m} \vec{F}$$
In pseudocode in 3D:
Timestep:
Let px, py, pz be the location of the object
at the start of the timestep
Let vx, vy, vz be the velocity of the object
at the start of the timestep
Let fx, fy, fz be the sum of forces acting on
the object during the timestep
Let m be the mass of the object
Let t be the duration of the timestep
Let c1 = 0.5*t*t/m
Let c2 = t/m
px = px + t*vx + c1*fx
py = py + t*vy + c1*fy
pz = pz + t*vz + c1*fz
vx = vx + c2*fx
vy = vy + c2*fy
vz = vz + c2*fz
# Timestep simulated for this object!
If the mass of the object changes during the timestep, for example like in a rocket or some kind of thruster, we should use
$$m(t) \vec{a}(t) = \vec{F}(t) + \vec{u}(t) \frac{d m}{d t}$$
where $u$ is the velocity vector of the exhaust, and $\frac{d m}{d t}$ is the rate of change of mass, negative for exhaust.
I do believe the continuous-time solution for $\vec{a}(t)$ (with constant force $\vec{F}$) is
$$\vec{a}(t) = \frac{\vec{F} + w \vec{u}}{m_0 + w t}$$
where mass changes as $m(t) = m_0 + w t$; that integrates to
$$\vec{v}(t) = \vec{v}(0) + \left ( \vec{u} + \frac{\vec{F}}{w} \right ) \log \left ( \frac{t w}{m_0} + 1 \right )$$
and again to
$$\vec{x}(t) = \vec{x}(0) + \left( \vec{F} \left ( \frac{t}{w} + \frac{m_0}{w^2} \right ) + \vec{u} \left ( t + \frac{m_0}{w} \right ) \right ) \log \left ( \frac{t w}{m_0} + 1 \right ) + t \left ( \vec{v}(0) - \vec{u} \right ) - \vec{F}\frac{t}{w}$$
At the limit $w = 0$, these simplify to the constant-mass functions earlier above. So, the pseudocode to implement a time step when the mass of the object changes, would be
Timestep:
Let px, py, pz be the location of the object
at the start of the timestep
Let vx, vy, vz be the velocity of the object
at the start of the timestep
Let fx, fy, fz be the sum of forces acting on
the object during the timestep
Let ux, uy, uz be the velocity of the exhaust
Let m be the mass of the object
at the start of the timestep
Let w be the rate of change of mass
(negative if exhaust)
Let t be the duration of the timestep
c0 = t / w
c1 = m / w
c2 = log(1 + t*w/m)
px = px + (fx*(c0 + c1/w + ux*(t + c1))*c2 + t*(vx - ux) + fx*c0
py = py + (fy*(c0 + c1/w + uy*(t + c1))*c2 + t*(vy - uy) + fy*c0
pz = pz + (fz*(c0 + c1/w + uz*(t + c1))*c2 + t*(vz - uz) + fz*c0
vx = vx + (ux + fx/w)*c2
vy = vy + (uy + fy/w)*c2
vz = vz + (uz + fz/w)*c2
m = m + t*w
# Timestep simulated for this object!