How to partition multivariate nonlinear function into a set of linear regions

Hi, I am studying control theory, specifically nonlinear models (using idnlgrey objects). The state update equation can be an arbitrary nonlinear function. A very simple example:

xnew = [x(1)^2+sin(x(2)); cos(x(1))^2+2*u];

This represents system behavior as a set of differential equations. There could be any number of states in the model.

I'd like to transform this continuous, nonlinear model into a set of continuous, linear models (approximating a nonlinear model as a set of linear regions). This requires breaking the function into regions, where within each region the state transformation is a linear mapping. Essentially redefine it in piecewise linear fashion, where each region has its own state-space equation, something like:

 xnew=A*x+B*u+e;

Does Matlab have a mechanism to support this? I know there is no "single" way to partition, but specifically I'm wanting hypercubes. And for now I don't have any requirement that the nonlinear->linear mapping is least-square error optimal....