The parameters enter nonlinearly into your equation. However, fortunately you can rewrite
your equation in such a way that you have linear parameters:
$$ (f(x+1) - f(x))^{1/3} = a f(x) - b $$
where $k_1 = a^3$ and $k_2 = b/a$.
So you can apply least squares to these linear equations in $a$ and $b$ for $x$ from $1$ to $n$.
EDIT: I tried some data that seem close to what you had in your picture:
$$ y = [237,130,120, 113, 111, 110] $$
Linear least squares for the residuals $(y_{i+1} - y_i)^{1/3} - a y_i + b$, $i = 1 \ldots,5$, in Maple produces
$$ a = -.0273525754168209,\ b = -1.67458695987193$$
which corresponds to
$$ k_1 = - 0.0000204641953284227606,\ k_2= 61.2222774036158200$$
I then tried nonlinear least squares for the residuals
$ y_{i+1} - y_i - k_1 (y_i - k_2)^3$, $i=1 \ldots 5$, using the above as initial values
(this is often necessary because the nonlinear least squares algorithms often provide only a local minimum rather than the global minimum). The result was
$$ k_1 =- 0.0000166852730695581124, k_2 = 51.1767746591434971 $$
The recursion using initial value $237$ and parameters from the linear least squares
produces the values
$$ 237, 125.855961638565, 120.330463855443, 116.104385196179, 112.721501867265, 109.926405798693$$
while the recursion using parameters from the nonlinear least squares produces
$$ 237, 129.938505804843, 121.786225982822, 115.912389815574, 111.385883230420, 107.744051206183$$
If you have some data $y$ and want to fit it to a curve of the form $1/(ax^2+bx+c)$, then what we can do is fit the data $1/y$ to the curve $ax^2+bx+c$. To do this in the least squares sense, we construct the following system of equations:
$$ \begin{pmatrix} x_1^2 & x_1 & 1 \\ x_2^2 & x_2 & 1 \\ \vdots & \vdots & \vdots \\ x_n^2 & x_n & 1\end{pmatrix} \begin{pmatrix} a \\ b \\ c\end{pmatrix} =
\begin{pmatrix}
1/y_1 \\ 1/y_2 \\ \vdots \\ 1/y_n
\end{pmatrix}
$$
The matix on the left-hand side is known as a Vandermonde matrix. Denote this matrix $A$
Since in general $n \neq 3$, to solve this we simply solve $A^TA \mathbf{p} = A^T\mathbf{y}$, where $\mathbf{p} = (a,b,c)^T$ is the parameter vector.
Thus, $$\mathbf{p} = (A^TA)^{-1}A^T\mathbf{y},$$ which may be computed in many ways (note that $A^TA$ is symmetric, so the LU factorization should be easy to compute).
EDIT: Regarding largeness of $\mathbf{y}$ (which should not matter for the data in your picture).
Note that we are solving $A^TA\mathbf{p} = A^T\mathbf{y}$. Look at $A^T\mathbf{y}$:
$$A^T \mathbf{y} = \begin{pmatrix}
\sum_{i=1}^n \frac{x_i^2}{y_1} \\
\sum_{i=1}^n \frac{x_i}{y_1} \\
\sum_{i=1}^n \frac{1}{y_1}
\end{pmatrix}$$
If $y_i$ is so small that you are concerned about numerical stability, solve the following problem instead:
$$10^m \hat{a} x^2 + 10^m \hat{b} x + 10^m \hat{c} = 1/y$$
for some $m > 0$. Then, you condition your Vandermonde matrix in such a way that $A^T\mathbf{y}$ becomes:
$$A^T \mathbf{y} = \begin{pmatrix}
\sum_{i=1}^n \frac{10^m x_i^2}{y_1} \\
\sum_{i=1}^n \frac{10^m x_i}{y_1} \\
\sum_{i=1}^n \frac{10^m}{y_1}
\end{pmatrix}$$
and the smallness is gone. Compute $\hat{a},\hat{b},\hat{c}$ in the traditional sense, then $a = 10^m\hat{a}$, and so forth.
Best Answer
Here's an example of what you're trying to do in Wolfram Alpha. I think you can do the same thing in Mathematica using similar, or identical notation.