Resolveré la ecuación de una forma genérica y así obtendré una fórmula.
Si la aceleración es
a = A - Bv^2 tendremos
$$\begin{align}&a=\frac{dv}{dt}=\frac{dv}{dx}·\frac{dx}{dt}=\frac{dv}{dx}v\\ & \\ & v \frac{dv}{dx}=A +Bv^2\\ & \\ & \frac {v}{A+Bv^2}dv = dx\\ & \\ & \int \frac {v}{A+Bv^2}dv = x+C\\ & \\ & \frac{1}{2B}\int \frac {2Bv}{A+Bv^2}dv = x+C\\ & \\ & \frac{ln(A+Bv^2)}{2B} = x+C\\ & \\ & ln(A+Bv^2) = 2B(x+C)\\ & \\ & A+Bv^2= e^{2B(x+C)}\\ & \\ & Bv^2= e^{2B(x+C)}-A\\ & \\ & v^2 = \frac{e^{2B(x+C)}-A}{B}\\ & \\ & v= \sqrt \frac{e^{2B(x+C)}-A}{B}\\ & \\ & \end{align}$$
Bueno, ya lo resolví con x en vez de y como en este caso. Además ya veré que tal me manejo con que la velocidad descendente sea positiva, para mí siempre ha sido negativa.
$$\begin{align}& \frac{ln(A+Bv^2)}{2B} = x+C\\ & \\ & \frac{ln(A+B·v_0^2)}{2B}=0+C\\ & \\ & C= \frac{ln(A+B·v_0^2)}{2B}\\ &\\ & v= \sqrt \frac{e^{2B\left(y+\frac{ln(A+B·v_0^2)}{2B}\right)}-A}{B}\\ &\\ &v=\sqrt{\frac{e^{2By}·e^{ln(A+Bv_0^2)}-A}{B}}\\ &\\ &v=\sqrt{\frac{e^{2By}(A+Bv_0^2)-A}{B}}\\ &\\ &\text{Y ssutituyendo tendremos}\\ &\\ &v(y)=\sqrt{\frac{e^{2·(-0.003)y}(9.81-0.003·3^2)-9.81}{-0.003}}\\ &\\ &v(y) = \sqrt{\frac{9.783e^{-0.006y}-9.91}{-0.003}}\\ &\\ &v(y)=\sqrt{\frac{9.91-9.783e^{-0.006y}}{0.003}}\\ &\\ &\text{Y la velocidad de régimen es cuando }\\ &\\ &y\to \infty\implies e^{-0.006y}\to 0\\ &\\ &v_r= \sqrt{\frac{9.91}{0.003}}\approx57.474632 m/s\end{align}$$
Y eso es todo.