I found this question very interesting since I had recently done some work on modeling projectile motion with drag.
Point 1: You are essentially updating the position and velocity using an explicit/forward Euler iteration where each new value for the states should be a function of the old values. In such a case, you should be updating the position first, then updating the velocity.
Point 2: There are more realistic physics models for the effect of drag friction. One model (suggested by Adam Liss) involves a drag force that is proportional to the velocity (known as Stokes' drag, which generally applies to low velocity situations). The one I previously suggested involves a drag force that is proportional to the square of the velocity (known as quadratic drag, which generally applies to high velocity situations). I'll address each one with regard to how you would deduce formulas for the maximum velocity and the time required to effectively reach the maximum velocity. I'll forego the complete derivations since they are rather involved.
Stokes' drag:
The equation for updating the velocity would be:
velocity += acceleration - friction*velocity
which represents the following differential equation:
dv/dt = a - f*v
Using the first entry in this integral table, we can find the solution (assuming v = 0 at t = 0):
v = (a/f) - (a/f)*exp(-f*t)
The maximum (i.e. terminal) velocity occurs when t >> 0, so that the second term in the equation is very close to zero and:
v_max = a/f
Regarding the time needed to reach the maximum velocity, note that the equation never truly reaches it, but instead asymptotes towards it. However, when the argument of the exponential equals -5, the velocity is around 98% of the maximum velocity, probably close enough to consider it equal. You can then approximate the time to maximum velocity as:
t_max = 5/f
You can then use these two equations to solve for f and a given a desired vmax and tmax.
Quadratic drag:
The equation for updating the velocity would be:
velocity += acceleration - friction*velocity*velocity
which represents the following differential equation:
dv/dt = a - f*v^2
Using the first entry in this integral table, we can find the solution (assuming v = 0 at t = 0):
v = sqrt(a/f)*(exp(2*sqrt(a*f)*t) - 1)/(exp(2*sqrt(a*f)*t) + 1)
The maximum (i.e. terminal) velocity occurs when t >> 0, so that the exponential terms are much greater than 1 and the equation approaches:
v_max = sqrt(a/f)
Regarding the time needed to reach the maximum velocity, note that the equation never truly reaches it, but instead asymptotes towards it. However, when the argument of the exponential equals 5, the velocity is around 99% of the maximum velocity, probably close enough to consider it equal. You can then approximate the time to maximum velocity as:
t_max = 2.5/sqrt(a*f)
which is also equivalent to:
t_max = 2.5/(f*v_max)
For a desired vmax and tmax, the second equation for tmax will tell you what f should be, and then you can plug that in to the equation for vmax to get the value for a.
This seems like a bit of overkill, but these are actually some of the simplest ways to model drag! Anyone who really wants to see the integration steps can shoot me an email and I'll send them to you. They are a bit too involved to type here.
Another Point: I didn't immediately realize this, but the updating of the velocity is not necessary anymore if you instead use the formulas I derived for v(t). If you are simply modeling acceleration from rest, and you are keeping track of the time since the acceleration began, the code would look something like:
position += velocity_function(timeSinceStart)
where "velocity_function" is one of the two formulas for v(t) and you would no longer need a velocity variable. In general, there is a trade-off here: calculating v(t) may be more computationally expensive than simply updating velocity with an iterative scheme (due to the exponential terms), but it is guaranteed to remain stable and bounded. Under certain conditions (like trying to get a very short tmax), the iteration can become unstable and blow-up, a common problem with the forward Euler method. However, maintaining limits on the variables (like 0 < f < 1), should prevent these instabilities.
In addition, if you're feeling somewhat masochistic, you may be able to integrate the formula for v(t) to get a closed form solution for p(t), thus foregoing the need for a Newton iteration altogether. I'll leave this for others to attempt. =)