state-space model of non-linear system?
A physical system is in state-space representation when we have a mathematical model of it as a set of input, output and state variables related by first-order differential equations only. 

The system 
$$ m\ddot{y} + b\dot{y} + k_1y+k_2y^3 = u$$ is not, since there's a second derivative. But, by introducing \$x_1 = \dot{y}\$, \$x_2 = y\$,
$$ \dot{\mathbf{x}} = \frac{d}{dt}\begin{pmatrix}\dot{y}\\y\end{pmatrix} =
\begin{bmatrix}\ddot{y}\\\dot{y}\end{bmatrix} = 
\left[\matrix{\frac{1}{m}(-b\dot{y}-(k_1 +k_2 y^2)y)\\\dot{y}}\right]+\left[\matrix{\frac{u}{m}\\0}\right] \\=
\left[\matrix{-\frac{b}{m}&amp;-\frac{k_1+k_2y^2}{m}\\1&amp;0}\right]\left[\matrix{\dot{y}\\y}\right]+\begin{bmatrix}\frac{1}{m}\\0\end{bmatrix} u,$$
or
$$ \dot{\mathbf{x}} = \mathbf{A}(\mathbf{x})\mathbf{x}+\mathbf{b}u $$
Notice that the system-matrix \$A\$ is not constant, but this is a state-space representation.

Update 1: Linearization

Given a non-linear system
$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}) + \mathbf{b}u,$$
the idea of linearization is to work out of a current set point \$\mathbf{x}_0,\$ so that
$$ \mathbf{f}(\mathbf{x}_0+\mathbf{d}) \approx \mathbf{f}(\mathbf{x}_0) + J_\mathbf{f}(\mathbf{x}_0)\cdot\mathbf{d}.$$

Here, we have 
$$\begin{eqnarray}f_1(x_1,x_2) &amp;=&amp; \frac{-bx_1 -k_1x_2 - k_2x_2^3}{m} \\ f_2(x_1,x_2) &amp;=&amp; x_1 \end{eqnarray}$$
so that the Jacobian matrix is
$$ J_f(\mathbf{x}) = \begin{bmatrix}\frac{\partial f_1}{\partial x_1} &amp; \frac{\partial f_1}{\partial x_2} \\ \frac{\partial f_2}{\partial x_1} &amp; \frac{\partial f_2}{\partial x_2}\end{bmatrix}=\begin{bmatrix}-\frac{b}{m} &amp; \frac{-k_1 -3k_2x_2^2}{m}\\ 1 &amp; 0\end{bmatrix}, $$

Notice that $$\dot{\mathbf{x}} = \frac{d}{dt} \mathbf{x}= \frac{d}{dt}(\mathbf{x}_0+\mathbf{d}) = \frac{d}{dt}\mathbf{d} = \dot{\mathbf{d}}.$$
So
$$ \dot{\mathbf{d}} = J_f(\mathbf{x}_0)\cdot\mathbf{d}+\mathbf{f}(\mathbf{x}_0)+\mathbf{b}u=\mathbf{Ad}+\mathbf{c}+\mathbf{b} u,$$
is a linear system that we can control normally, with
$$ \mathbf{A} = J_F(\mathbf{x}_0),$$ and 
$$ \mathbf{c} = \mathbf{f}(\mathbf{x}_0). $$

Update 2: Lyapunov Control

It looks like your problem is similar to a non-linear spring model. For instance, a system with a hardering spring is often modeled 
$$ F_s(x) = k(1+a^2y^2)y$$
which corresponds your equation with \$k_1 = k, k_2 = ka^2\$. A way you can solve non-linear control is by Lyapunov theory. For example, let us assume that we have a desired input signal \$y_d\$, and that we want to ensure that the error $$e = y_d-y \to 0 $$
Define the function $$ r = \dot{e}+\alpha e,$$ which gives a weighted notion of error of both position and velocity, and define the "energy" function $$ V = \frac{1}{2} r^2,$$ which is strictly positive for all \$\dot{y} \ne 0\$ and \$y\ne 0\$. Look at the time derivative
$$ \dot{V} = r\;\dot{r}= (\dot{e}-\alpha e)(\ddot{e}+\alpha\dot{e})= (\dot{e}-\alpha e)(\ddot{y_d}-\ddot{y}+\alpha\dot{e})$$
The goal is to make sure that the energy is exponentially decreasing, thus 
$$ \dot{V} = -\kappa V,$$ so we want 
$$ (\dot{e}-\alpha e)(\ddot{y}_d - \ddot{y} + \alpha \dot{e}) = -\kappa \frac{1}{2}(\dot{e}-\alpha e)^2, $$ or
$$ (\ddot{y}_d - \ddot{y} + \alpha \dot{e}) = -\frac{\kappa}{2}(\dot{e}-\alpha e). $$
Now we can insert the system dynamic for \$\ddot{y}\$ to get
$$ \ddot{y}_d - \frac{u-b\dot{y}-k_1y-k_2y^3}{m} + \alpha \dot{e} = -\frac{\kappa}{2}(\dot{e}-\alpha e),$$
and solving for \$u\$, we get the control equation:
$$ u-b\dot{y}-k_1y-k_2y^3 = m\left(\alpha \dot{e} + \ddot{y}_d + \frac{\kappa}{2}(\dot{e}-\alpha e)\right)\\ u = b\dot{y}+k_1y+k_2y^3 + m\alpha \dot{e} + m\ddot{y}_d + \frac{m\kappa}{2}\dot{e}-\frac{m\kappa\alpha}{2}e, $$
with tunable parameters.