Two principles here:
- When dealing with a differential equation, you define intermediate state variables so everything is in terms of first derivatives.
- This system is nonlinear, so the state-space equations won't be in terms of matrices.
Applying these principles, we define a state vector:
$$
\mathbf x = [x_1, x_2]^T,
$$
where:
$$
x_1 = y \\
x_2 = \dot y
$$
Note that:
$$
\dot x_1 = x_2
$$
Substituting these into your original equation yields:
$$
\dot x_1 = x_2 = f_1(x_2) \\
\dot x_2 = g - \frac{C}{m} \frac{u^2}{x_1^2} = f_2(x_1, u)
$$
which is of the form:
$$
\dot{\mathbf{x}} =
\mathbf f(\mathbf x, u, t)
$$
Since the model coefficients (i.e. $C, g, m$) don't depend on time, we can drop the $t$:
$$
\dot{\mathbf{x}} =
\mathbf f(\mathbf x, u)
$$
Now that you have a put the ODE in explicit form, you can take derivatives to find the operating point.
The basic idea is to approximate $\mathbf f(\mathbf x, u)$ with a first-order Taylor approximation about some operating point $(\mathbf x_0, u_0)$:
$$
\mathbf f(\mathbf x, u) \approx \mathbf f(\mathbf x_0, u_0) + (\mathbf x - \mathbf x_0) \frac{\partial \mathbf f}{\partial \mathbf x}\Bigg|_{(\mathbf x_0, u_0)} + (u - u_0) \frac{\partial \mathbf f}{\partial u} \Bigg|_{(\mathbf x_0, u_0)}
$$
EDIT:
From the referenced lecture slides, note that:
$$
\mathbf f(\mathbf x, u) = \mathbf f(\mathbf x_0 + \Delta \mathbf x, u_0 + \Delta u) = \dot{\mathbf{x}}|_{(\mathbf x_0, u_0)} + \dot{\Delta \mathbf{x}},
$$
where $\Delta \mathbf x = \mathbf x - \mathbf x_0$, $\Delta u = u - u_0$, and
$$
\dot{\mathbf{x}}|_{(\mathbf x_0, u_0)} = f(\mathbf x_0, u_0).
$$
This leads to a linear, time-invariant system in terms of the state variable $\Delta \mathbf x$:
$$
\dot{\Delta \mathbf{x}} = \mathrm A \Delta \mathbf x + \mathrm B \Delta u,
$$
where:
$$
\mathrm A = \frac{\partial \mathbf f}{\partial \mathbf x}\Bigg|_{(\mathbf x_0, u_0)} \\
\mathrm B = \frac{\partial \mathbf f}{\partial u}\Bigg|_{(\mathbf x_0, u_0)}
$$
The solution for this coupled set of equations yields the displacements $\Delta \mathbf x$ from the operating point $(\mathbf x_0, u_0)$. This is true for any operating point $(\mathbf x_0, u_0)$, regardless of whether or not $\mathbf f(\mathbf x_0, u_0) = 0$. Also note that this system is LTI mathematically, not physically. In reality, the matrices $\mathrm A$ and $\mathrm B$ can also change with time.
Note that a good control algorithm based on a linearized model will have to update the operating point at a fast enough rate that the nonlinearities don't cause any problems.