I'm going to give two answers.
Answer 1: If you already know $a$, then $c = \bar{y} - a\bar{x}$, where $\bar{x}$ and $\bar{y}$ are the means of the $x$ and $y$ values, respectively. (This is actually in that pile of formulas in the link provided by Yuval Filmus.)
Answer 2: At least in my version of Matlab, polyfit gives both $a$ and $c$. For example, if you generate a set of data via (I'm modifying the example given in the Matlab help)
x = (0: 0.1: 2.5)';
y = erf(x);
and then call polyfit with $n=1$ (for a linear curve fit, or linear regression)
p = polyfit(x,y,1)
The output is
0.3554 0.3191
which means your linear equation is $y = 0.3554x + 0.3191$.
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$$
Best Answer
Problem statement
Start with a polynomial and its derivatives: $$ \begin{align} y(x) &= a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3}, \\ y'(x) &= a_{1} x + 2 a_{2} x + 3 a_{3} x^{2}, \\ y''(x) &= 2 a_{2} + 6 a_{3} x. \end{align} $$ (Notice that we are an affine transformation from any other set of polynomial functions, such as those of Legendre.) The fit is looking for the best set of $n=4$ parameters $a$.
There are a sequence of measurements: one set measures functions values, another derivatives, the last, the second derivative. $$ \left\{ x_{1,k}, y_{k} \right\}_{k=1}^{\mu_{1}}, \quad \left\{ x_{2,k}, y'_{k} \right\}_{k=1}^{\mu_{2}}, \quad \left\{ x_{3,k}, y''_{k} \right\}_{k=1}^{\mu_{3}}. $$ This general method accounts for different types of measurements (function value, first derivative, second derivative) at different locations $x$.
Construct linear system
Your conditions lead to a block structure for the system matrix $\mathbf{A}$. The problems have different measurement locations $x$ and measurements $y$, but they share the amplitudes. $$ \begin{align} \mathbf{A} a &= y \\ % \left[ \begin{array}{cccc} % 1 & x_{1,1} & x_{1,1}^{2} & x_{1,1}^{3} \\ 1 & x_{1,2} & x_{1,2}^{2} & x_{1,2}^{3} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1,\mu_{1}} & x_{1,\mu_{1}}^{2} & x_{1,\mu_{1}}^{3} \\\hline % 0 & 1 & 2 x_{2,1} & 3x_{2,1}^{2} \\ 0 & 1 & 2 x_{2,2} & 3x_{2,2}^{2} \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 1 & 2 x_{2,\mu_{2}} & 3x_{2,\mu_{2}}^{2} \\\hline % 0 & 0 & 2 & 6x_{3,1} \\ 0 & 0 & 2 & 6x_{3,2} \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 2 & 6x_{3,\mu_{3}} % \end{array} \right] % \left[ \begin{array}{c} a_{0} \\ a_{1} \\ a_{2} \\ a_{3} \end{array} \right] % & = % \left[ \begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{\mu_{1}} \\\hline y'_{1} \\ y'_{2} \\ \vdots \\ y'_{\mu_{2}} \\\hline y''_{1} \\ y''_{2} \\ \vdots \\ y''_{\mu_{3}} \end{array} \right] % \end{align} $$
The data vector has $m= \mu_{1} + \mu_{2} + \mu_{3}$ rows, and the solution vector is of length $n=4$. In general, the problem will have full column rank $\rho = n$. The dimensions are $$ \mathbf{A} \in \mathbb{C}^{m \times n}_{\rho}, \quad y \in \mathbb{C}^{n}, \quad a \in \mathbb{C}^{m} $$
Least squares solution
The least squares minimizers are the defined as $$ a_{LS} = \left\{ a \in \mathbb{C}^{n} \colon \lVert \mathbf{A} a - y \rVert_{2}^{2} \text{ is minimized} \right\} $$ A rich toolkit offers many paths to solution. One method is the normal equations $$ \mathbf{A}^{*} \mathbf{A} a = \mathbf{A}^{*} y $$ which has the solution $$ a = \left( \mathbf{A}^{*} \mathbf{A} \right)^{-1} \mathbf{A}^{*} y. $$
Example
Start with ideal data and add random noise. The solution vector $$ % a = % \left[ \begin{array}{r} a_{0} \\ a_{1} \\ a_{2} \\ a_{3} \end{array} \right] % = % \left[ \begin{array}{r} 1 \\ -2 \\ 3 \\ 4 \end{array} \right]. % $$ The data sets are $$ % x_{1} = % \frac{1}{10} \left[ \begin{array}{r} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \\ 7 \\ 8 \\ 9 \\ 10 \end{array} \right], \, % y = % \frac{1}{500} \left[ \begin{array}{r} 417 \\ 376 \\ 389 \\ 468 \\ 625 \\ 872 \\ 1221 \\ 1684 \\ 2273 \\ 3000 \end{array} \right], \quad %%% x_{2} = % \frac{1}{6} \left[ \begin{array}{r} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \end{array} \right], \, % y' = % \frac{1}{3} \left[ \begin{array}{r} -2 \\ 4 \\ 12 \\ 22 \\ 34 \\ 48 \end{array} \right], \quad %%% x_{3} = % \frac{1}{5} \left[ \begin{array}{r} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{array} \right], \, % y'' = % \left[ \begin{array}{r} 10 \\ 14 \\ 18 \\ 22 \\ 26 \\ 30 \end{array} \right] $$ Before perturbation, the linear system looks like this: $$ \left[ \begin{array}{rrrr} 1 & 0.1 & 0.01 & 0.001 \\ 1 & 0.2 & 0.04 & 0.008 \\ 1 & 0.3 & 0.09 & 0.027 \\ 1 & 0.4 & 0.16 & 0.064 \\ 1 & 0.5 & 0.25 & 0.125 \\ 1 & 0.6 & 0.36 & 0.216 \\ 1 & 0.7 & 0.49 & 0.343 \\ 1 & 0.8 & 0.64 & 0.512 \\ 1 & 0.9 & 0.81 & 0.729 \\ 1 & 1 & 1 & 1. \\\hline 0 & 1 & 0.333333 & 0.0833333 \\ 0 & 1 & 0.666667 & 0.333333 \\ 0 & 1 & 1 & 0.75 \\ 0 & 1 & 1.33333 & 1.33333 \\ 0 & 1 & 1.66667 & 2.08333 \\ 0 & 1 & 2 & 3 \\\hline 0 & 0 & 2 & 1.2 \\ 0 & 0 & 2 & 2.4 \\ 0 & 0 & 2 & 3.6 \\ 0 & 0 & 2 & 4.8 \\ 0 & 0 & 3 & 6 \end{array} \right] % \left[ \begin{array}{r} a_{0} \\ a_{1} \\ a_{2} \\ a_{3} \end{array} \right] % = % \frac{1}{1500} \left[ \begin{array}{r} 1251 \\ 1128 \\ 1167 \\ 1404 \\ 1875 \\ 2616 \\ 3663 \\ 5052 \\ 6819 \\ 9000 \\\hline -1000 \\ 2000 \\ 6000 \\ 11000 \\ 17000 \\ 24000 \\\hline 16200 \\ 23400 \\30600 \\ 37800 \\ 45000 \end{array} \right] % $$ A random error $-0.1 < \epsilon < 0.1$ was added to each $y$ value.
The normal equations are $$ % \frac{1}{1800000} \left[ \begin{array}{rrrr} 18000000 & 9900000 & 6930000 & 5445000 \\ 9900000 & 17730000 & 18045000 & 18209940 \\ 6930000 & 18045000 & 58759940 & 90824850 \\ 5445000 & 18209940 & 90824850 & 174558629 \\ \end{array} \right] % a = % \left[ \begin{array}{r} 23.2848 \\ 56.2244 \\ 283.846 \\ 522.72 \end{array} \right] % $$ The perturbed solution is $$ a_{LS} = \left( \mathbf{A}^{*} \mathbf{A} \right)^{-1} \mathbf{A}^{*} y = \left[ \begin{array}{r} 1.10766 \\ -2.09876 \\ 3.02645 \\ 3.99983 \end{array} \right] $$