You have a system $\dot x=C(k)T(x)$ where $C(k)$ is a matrix containing the parameters $k$ and $T(x)$ are terms in $x$ without parameters.
If you were able to measure sufficiently fine time series of all the components to then replace the derivatives $\dot x(t_k)$ as difference quotients (or the slope of local linear regression) $v_k$, then the resulting system $v_k=C(k)T(x_k)$ is indeed linear in the parameters.
However, you most likely have only one or two components measured. The parameters in the equation for the other components do only influence the measured components over the non-linear dependence. Thus it is necessary to solve the system numerically for a current approximation of the parameters, compute derivatives/sensitivities, etc. In the time propagation of the solver, the parameters are also propagated through the non-linear dynamics, thus will occur non-linearly in the predicted component values. Even in the easiest case of the scalar equation $\dot x=kx$, the solution $x=Ce^{kt}$ is non-linear. Thus a non-linear problem.
Or you could look at the Euler method as most simple example. While $x_1=x_0+hC(k)T(x_0)$ is still linear,
$$
x_2=x_0+hC(k)T(x_0)+hC(k)T(x_0+hC(k)T(x_0))
$$
already contains the parameters in quadratic or higher order. One can observe that $hC(k)$ is always combined, so that a degree $d$ term in $k$ has the coefficient $h^d$, but the number of terms of a given degree will also grow superlinearly.
As sort of a corollary, the shorter the integration time span, the less important will the non-linear terms in the parameters be. Splitting a fitting task into smaller sub-tasks with such nearly linear dependence on the coefficiens can help to improve and stabilize the convergence to useful parameters.
