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
