In fact, Newton's law works in terms of $$\vec F=\frac {d\vec p}{dt}$$ where $\vec F$ is the resulting external force and $\vec p$ is momentum.
The expression for momentum writes $m\vec v$, where $m$ is the mass and $\vec v$ is the speed. In most cases the mass is constant, hence we simplify $\frac {d\vec p}{dt}=m\frac {d\vec v}{dt}$.
However, in your case mass is not constant, so we are forced to use the equation $$\vec F=m(t)\frac {d\vec v(t)}{dt}+\frac {dm(t)}{dt}\vec v.$$
In order to graph the projectile motion of a sphere, we must be able to calculate the forces acting on said sphere at a given time. One of these forces is drag, $F_d(t)$. $$F_d(t) = \dfrac{1}{2} \cdot C_d \cdot \rho \cdot v(t) \cdot |v(t)| \cdot A$$
Where $C_d$ is the drag coefficient, $\rho$ is the density of air, $v(t)$ is the velocity vector at time $t$ and $A$ is the cross sectional area of the sphere. Drag is a two dimensional vector such that $$F_d(t) = \begin{bmatrix} F_d^x (t) \\ F_d^z (t)\end{bmatrix}$$
We can therefore work out the sum of forces acting in $x$ and $z$ by using $\sum f = ma$. In this case, $\sum f^x (t) = -F_d^x (t) = ma^x (t)$. As such, $$a^x (t) = \dfrac{-F_d^x (t)}{m}$$
When considering the $z$ component, $\sum f^z (t) = -F_d^z (t) - mg = ma^z (t)$. Therefore, $$a^z(t) = \dfrac{-F_d^z (t)}{m} - g$$
According to Euler's method, we can add a small amount of time $\delta$ (say $0.0005$) and calculate the coordinates at $t + \delta$ where $t$ is the current time. Therefore,
$$v^x (t + \delta) = v^x (t) + a^x (t) \cdot \delta$$
$$v^z (t + \delta) = v^z (t) + a^z (t) \cdot \delta$$
We can then use these formulae to work out the coordinate (displacement) vector, $s(t)$ at time $t$.
$$s^x (t + \delta) = s^x (t) + v^x (t) \cdot \delta$$
$$s^z (t + \delta) = s^z (t) + v^z (t) \cdot \delta$$
We now almost have all of the information we need to begin calculating the coordinates of the motion of the sphere. Although, we must first define some values. $C_d = 0.5$ as this is the standard drag coefficient for a sphere, $\rho = 1.225$ as this is the density of air, let the mass $m = 1.0$ and $A = \pi r^2$ where $r = 0.5$. We know that $c(0) = \begin{bmatrix} 0 \\ 0\end{bmatrix}$, and we can define the initial velocities as $u(0) = \begin{bmatrix} 50 \\ 30\end{bmatrix}$.
We can therefore calculate $F_d(0)$ which will allow us to calculate $a^x (0)$, $a^z(0)$ meaning we can calculate $v(0 + \delta)$, and therefore $c(0 + \delta)$. This process should be repeated until $s^z(t) \leq 0$, as this means we have reached or gone below the $z$ axis. Once this has occurred, we can plot our coordinates and set a range of $z \geq 0$, and our graph will be complete.
Note that we have ignored the $y$ axis as the sphere is being projected directly along the $x$ axis, and so changes only occur in $x$ and $z$.
Thanks to @David K for his help.
Below is a proof for the Euler's method equations.
Consider the acceleration vector $a(t)$. Since $a(t) = \dfrac{dv}{dt}$, $a(t) \approx \dfrac{\Delta v}{\Delta t}$ = $\dfrac{v(t + \delta) - v(t)}{\delta}$. Therefore, $$v(t + \delta) = v(t) + a(t) \cdot \delta$$
We can then use this formula to work out the displacement vector, $s(t)$ at time $t + \delta$ which will give us our new coordinates. Consider the velocity vector $v(t)$. As $v(t) = \dfrac{ds}{dt}$, $v(t) \approx \dfrac{\Delta s}{\Delta t} = \dfrac{s(t + \delta) - s(t)}{\delta}$. Therefore, $$s(t + \delta) = s(t) + v(t) \cdot \delta$$
Best Answer
Substitute $m=\frac{4\pi\delta}3r^3$ into the Newton’s law $F=mg = \frac{d}{dt}(mv)$ to get $$r^3g = \frac{d}{dt}(r^3v)\tag1$$ and also into the proportional condensation $\frac{dm}{dt}= 4\pi r^2 c$ to get $$\frac{dr}{dt}=\frac c\delta$$ Then, plug the resulting solution $r= \frac c\delta t$ into (1) to get $$t^3g =\frac{d}{dt}(t^3v)$$ Integrate both sides to obtain $\frac14 t^4g = t^3 v$, or $v=\frac14gt$, i.e. accelerating at $\frac14g$.