I'm not sure if there even exists a closed form for drag or wind, but it is quite easy to simulate in a step-wise fashion (like all the physics libraries do):
set your initial condition:
$$ x, y, v_x, v_y \; (\text{for }t=0) $$
update position:
$$ x = x + (v_x \times dt) \\ y = x + (v_y \times dt) $$
(where dt is the time elapsed since the last update, aka delta time)
calculate these velocity helpers:
$$ \begin{align} v^2 &= (v_x)^2 + (v_y)^2 \\ \lvert v \rvert &= \sqrt{v^2} \end{align} $$
(where \$\lvert v \rvert\$ represents the length of \$v\$)
calculate drag force:
$$ f_{drag} = c \times v^2 $$
(where c is the coefficient of friction small!)
accumulate forces:
$$ \begin{align} f_x &= \left(-f_{drag} \times {v_x \over \lvert v \rvert}\right) \\ f_y &= \left(-f_{drag} \times {v_y \over \lvert v \rvert}\right) + (-g \times mass) \end{align} $$
(where \$mass\$ is the mass of your golf ball)
update velocity:
$$ v_x = v_x + f_x \times \frac{dt}{mass} \\ v_y = v_y + f_y \times \frac{dt}{mass} $$
That's basically Euler's Method for approximating those physics.
A bit more on how the simulation as requested in the comments:
- The initial condition \$(t = 0)\$ in your case is
$$ \begin{align} x &= 0 \\ y &= 0 \\ v_x &= v_0 \times cos(\theta) \\ v_y &= v_0 \times sin(\theta) \end{align} $$
It's basically the same as in your basic trajectory formula where every occurrence of t is replaced by 0.
The kinetic energy \$KE = 0.5m(V^2) \$ is valid for every \$t\$. See \$v^2\$ as in (3) above.
The potential energy \$ PE = m \times g \times y \$ is also always valid.
If you want to get the current \$(x,y)\$ for a given \$t_1\$, what you need to do is initialize the simulation for \$t = 0\$ and do small dt updates until \$t = t_1\$
If you already calculated \$(x,y)\$ for a \$t_1\$ and you want to know their values for a \$t_2\$ where \$t_1 \lt t_2\$, all you need to do is calculating those small dt update steps from \$t_1\$ to \$t_2\$
Pseudo-Code:
simulate(v0, theta, t1) dt = 0.1 x = 0 y = 0 vx = v0 * cos(theta) vy = v0 * sin(theta) for (t = 0; t < t1; t += dt) x += vx * dt y += vy * dt v_squared = vx * vx + vy * vy v_length = sqrt(v_squared) f_drag = c * v_squared f_grav = g * mass f_x = (-f_drag * vx / v_length) f_y = (-f_drag * vy / v_length) + (-f_grav) v_x += f_x * dt / mass v_y += f_y * dt / mass end for return x, y end simulate
