Linear system with non-constant matrix containg entries form unknown vector

Consider a system of PDEs $$\begin{cases} u_t = \nabla \cdot (D(u)\nabla u) + \frac{c}{K_U+c}u-ku\\ c_t = d_c\Delta c -\frac{\nu_U c}{K_U + c}u \end{cases}$$ with some boundary conditions. Here, $$D(u)$$ is a diffusion coefficient which depends on $$u$$; $$K_U$$, $$\nu_U$$ and $$k$$ are some constants. $$D(u)$$ can be defined as, for example,

$$D(u):=\delta \frac{u^\alpha}{(1-u)^\beta},$$ with $$\alpha,\beta,\delta$$ being some constants.

After this system of PDEs is discretized using the finite volume method one obtains the following system of ODEs

$$\begin{cases}\frac{d\vec{U}}{dt}=\underline{D(\vec{U})}\vec{U}+\underline{R_U(\vec{C})}\vec{U}\\ \frac{d\vec{C}}{dt}=\underline{L}\vec{C}-\underline{R_C(\vec{C})}\vec{U}+\vec{b} \end{cases}$$

where the underlined letters are matrices and $$\vec{b}$$ is a vector containing some terms from (unspecified here) the boundary conditions.

As we can see, the matrix $$\underline{D(\vec{U})}$$ depends on the vector for which a numerical method will solve this system of equations, that is the vector $$\vec{U}$$.

But then how can such a system be solved linearly if it will contain non-linear terms? I.e. what am I missing?