Linearize differential equation Given the following differential equation:
$$
\dot{y}_1 + p(t)f(y_1) = g(t)
$$
then I want to linearize it in the sense that it becomes a system of equations on the form:
$$
\dot{y}_1 + p_1(t)y_1(t) = g_1(t)
$$
$$
\dot{y}_2 + p_2(t)y_2(t) = g_2(t)
$$
$$ \vdots $$
$$
\dot{y}_n + p_n(t)y_n(t) = g_n(t)
$$
where I allow for any finite number of additional linear differential equations, also the functions $g_i(t)$ and $p_i(t)$ may depend on any $y_j$ where $i\neq j$

Is this possible? preferably exact (of course) or at least as an approximation where it converges towards the exact case for $n\to \infty$ Any ideas or references are highly appreciated.
No, that is in general not possible. It would be a miracle making decades of research into numerical methods worthless. 

What you can do is linearize along a solution, so that a close solution $\bar y=y_0+εy_1$ satisfies in the first order the linearized system
\begin{align}
\dot y_0+f(y_0)&amp;=g\\
\dot y_1+f'(y_0)·y_1&amp;=0\\
\end{align}
One might expand to higher order terms by including higher orders in the perturbation parameter $ε$ in $\bar y=y_0+εy_1+ε^2y_2+…$ which then includes non-linear terms in the lower order coefficients.
Warning: this answer may contain traces of engineering language and therefore lack of rigor.

One way to do this would be to split the functions into representations which are valid at different intervals where we assume the function to be nice in each local neighbourhood. Then we require the function values at the boundaries between these sets to be smooth. We could for instance be using polynomials windowed with a function from the exponential family.

Easiest example would probably be a gaussian which windows a set of polynomials. They are easy to differentiate so we can easily get linear constraints on the polynomial coefficients.