So finally, this should yield $\Delta v = \sqrt{\frac{2GM}{r_o}} - \sqrt{\frac{2GM}{r_f}}$ plus any extra speed we wanted to have at that height, maybe enough to sustain a stable orbit. Is my reasoning correct?
No, it's not. You are calculating the difference between escape velocity at two different altitudes.
To put something into orbit, the object needs to gain kinetic and potential energy. The amount of energy expended (delta V) depends very much on the path taken to achieve this goal. Suppose we could teleport to the desired orbit, a circular orbit some height $h$ above the surface of an airless, non-rotating planet. The teleport device of course obeys conservation of energy. The initial specific energy is $E_0 = -\frac \mu r$, where $\mu \equiv GM$ is the planet's standard gravitational parameter. The final specific energy is $E_1 = \frac {v^2} 2 - \frac \mu {r+h} = -\frac 1 2 \frac \mu{r+h}$. The difference is the minimum specific energy needed to achieve that orbit:
$$\Delta E = \frac\mu r-\frac 12\frac\mu {r+h} = \frac\mu2\left(\frac2r - \frac1{r+h}\right)$$
We can't teleport, of course. Any real scheme will require more energy than this. Not quite real yet, suppose the vehicle applies an impulsive (instantaneous) delta V at the surface of the planet, with the velocity vector angled at some angle $\phi$ above the horizontal and the magnitude just enough to have the object reach a height $h$ above the surface at aphelion. Then the vehicle applies another impulsive burn to circularize the orbit.
The semi-major axis, initial velocity, aphelion velocity, and eccentricity of the transfer orbit (the path from the surface to the aphelion point) are
$$\begin{aligned}
a &= \frac{r+h}{1+\epsilon} \\
{v_i}^2 &= \mu \left(\frac2r-\frac1a\right) \\
{v_a}^2 &= \mu \frac{1-\epsilon}{r+h} \\
\epsilon &= \frac{h^2 + r(r+2h)\sin^2\phi}{h(2r+h)+r^2\sin^2\phi}
\end{aligned}$$
The first is simple geometry. The aphelion distance is $a(1+\epsilon)$, and we want that to be $r+h$. The second and third are the vis-viva equation, which comes in very handy for these kinds of calculations. The final expression is a tedious calculation involving the eccentricity vector, $\vec \epsilon = \frac {\vec v \times (\vec r \times \vec v)}\mu - \hat r$, the semi-major axis, and the initial velocity.
An example using $\mu = 398600.4418\,\text{km}^3/\text{s}^2$ (the Earth's gravitational parameter), $r = 6378.1366\,\text{km}$ (the Earth's equatorial radius), $h=2000\,\text{km}$ (the altitude in the question), and various angles follows.
phi(°) ecc a(km) v1(km/s) v2(km/s) dv2(km/s) dvtot(km/s) loss(km/s)
0 0.135536 7378.137 8.424077 6.413111 0.484444 8.908522 0.110036
30 0.517815 5519.867 7.264824 4.789631 2.107924 9.372747 0.574261
60 0.893738 4424.125 5.906997 2.248450 4.649105 10.556103 1.757617
90 1.0 4189.068 5.462334 0.0 6.897555 12.359889 3.561403
In the above, phi is the angle above horizontal, ecc is the eccentricity of the transfer orbit, a is the semi-major axis of the transfer orbit, v1 is the initial velocity needed to achieve that transfer orbit, v2 is the velocity on the transfer orbit at aphelion (which is occurs at 2000 km above the surface of the Earth), dv2 is the delta V needed to circularize the orbit at aphelion, dvtot is the sum of v1 and dv2, and loss is amount by which dvtot exceeds the ideal teleport delta V (8.798486 km/s, in this example).
A key takeaway point: On an airless planet (or moon), it makes much more sense to launch horizontally than vertically. Doing so from the surface of the Earth makes no sense. This would mean flying at Mach 25 through the thickest part of the atmosphere. Rockets launch vertically because the first goal is to get above the thickest part of the atmosphere. Rockets then start a turn toward the horizontal so as to avoid an overly huge penalty for flying vertically.
I used a very simple python script to generate the above table. Python is a handy language for such analyses. Here's the script:
import math
mu = 398600.4418 # Earth gravitational parameter, in km^3/s^2
r = 6378.1366 # Earth equatorial radius, in km (IERS 2003)
h = 2000.0 # Desired altitude, in km
vc = math.sqrt(mu/(r+h))
dv_opt = math.sqrt(mu*(2.0/r - 1.0/(r+h)))
print(dv_opt)
for f in range(0,91,30):
theta = f * math.pi/180.0
sintsq = math.sin(theta)**2
ecc = (h*h+(r*r+2*r*h)*sintsq)/(h*h+2*r*h+r*r*sintsq)
a = (r+h)/(1.0+ecc)
v1 = math.sqrt(mu*(2.0/r-1.0/a))
v2 = math.sqrt(mu*(2.0/(r+h)-1.0/a))
dv = v1 + (vc - v2)
print("{:2d} {:.6f} {:8.6f} {:8.6f} {:8.6f} {:8.6f} {:9.6f} {:8.6f}"
.format(f, ecc, a, v1, v2, vc-v2, dv, dv-dv_opt))
