First order expansion of Euler-Lagrange equations

I know that in field theory Euler Lagrange equations are $$p_i-d_\mu p^\mu_i=0$$. (Classical notations, $$p_i=\frac{\partial L}{\partial y^i}, p_i^\mu=\frac{\partial L}{\partial y^i_\mu}$$). Being a differential equation, I can expand it to first (or higher) order, but how can we do this?

I was wondering if it is right to consider the expansion of a solution as $$y^i=y_0^i+\epsilon y^i_1+...$$ and then substitute this expression in the E.L. equations to get the first order expansion?

I tried to do the computation but I wasn't able to get any result, any suggestion?

I don't recommend to use the notation $$d_\mu$$ since the derivative that appears in the E.L equations of motion is partial derivative $$\partial_\mu$$. Then if the Lagrangian density depends on the field $$y$$ and its derivatives $$\partial_\mu y$$ the E.L equations of motion read
$$$$\partial_\mu\dfrac{\partial \mathcal L}{\partial (\partial_\mu y)}-\dfrac{\partial \mathcal{L}}{\partial y}$$=0$$
$$y=c_0y_0+c_1 y_1+c_2y_2+\ldots=\displaystyle\sum_{i=1}^{\infty}c_i y_i$$
Then note that the Lagrangian density becomes a function of that many fields and its derivatives $$\mathcal{L}(\partial_\mu y, y)\equiv \mathcal{L}(\partial_\mu y_0,\partial_\mu y_1,\ldots, y_0,y_1,\ldots)$$. So that when you derive (through the variational principle) the Euler-Lagrange equations of motion you get a system of (possibly coupled) differential equations for the $$y_0,y_1,\ldots$$ fields.
$$$$\displaystyle\sum_{i=1}^{\infty}\bigg[\partial_\mu\dfrac{\partial \mathcal L}{\partial (\partial_\mu y_i)}-\dfrac{\partial \mathcal{L}}{\partial y_i}\bigg]$$=0.$$