Let us consider the following nonlinear polynomial system $$\dot{x} = f(x,u),$$ where $x=[x_1, ... , x_n]$. A Taylor expansion about $(x_0,u_0)$ gives
$$f(x,u) = f(x_0,u_0) + \frac{\partial f}{\partial x}|_{x=x_0,u=u_0}\cdot (x-x_0) + \frac{\partial f}{\partial u}|_{x=x_0,u=u_0}\cdot (u-u_0)+r(x,u)$$
where $r(x,u)$ is the remainder.
Here, $f(x_0,u_0)$ can be found by steady state conditions $\dot{x}=0=f(x_0,u_0)$
However, the result of linearizing with a high order >2 is that $x$ can still be present in a nonlinear way. For example: $\frac{\partial (x_1x_3^4)}{\partial x_1}|_{x_1=x_{1,0},u=u_0} = 4x_3^3$
I am confused how to tackle this, as the function is still nonlinear, while I have linearized. Could someone please explain how to do this?
Welcome to MSE!
I think you confused two things, the partial derivative as a function of $x$ and the partial derivative as numbers / matrix (which is the evaluation of function at a point)
In your equation $$ f(x,u) = f(x_0,u_0) + \frac{\partial f}{\partial x}|_{x=x_0,u=u_0}\cdot (x-x_0) + \frac{\partial f}{\partial u}|_{x=x_0,u=u_0}\cdot (u-u_0)+r(x,u) $$ $\frac{\partial f}{\partial x}|_{x=x_0,u=u_0}$ is a number / matrix that you evaluated at $x=x_0,u=u_0$. This is not a function of $x$
Thus, we can say the approximation $$ f(x_0,u_0) + \frac{\partial f}{\partial x}|_{x=x_0,u=u_0}\cdot (x-x_0) + \frac{\partial f}{\partial u}|_{x=x_0,u=u_0}\cdot (u-u_0) $$ is a linear function in $x$ and $u$, thus it linearized the original $f(x,u)$