Is propagation of uncertainties linear?

It's not, at least not in the statistical sense. What you are doing is finding the (linearly approximated) change in $y$ obtained by changing the inputs by their standard deviations. This is okay as a rough approximation, but you can do better with almost no extra work.

If you want the actual variance and standard deviation of $y$, the formula is different. Suppose $$ y = f(x_1, \ldots, x_n) = f_0 + \sum_{i=1}^n a_i x_i + \sum_{i=1}^n \sum_{j=i}^n b_{ij} x_i x_j + \mathcal{O}(x^3). $$ First, we can compute the expectation of $y$, $E(y) \equiv \hat{y}$, very easily: $$ \hat{y} = f_0 + \sum_{i=1}^n a_i \hat{x}_i + 2 \sum_{i=1}^n \sum_{j=i+1}^n b_{ij} \hat{x}_i \hat{x}_j + \sum_{i=1}^n b_{ii} E(x_i^2) + \mathcal{O}(x^3). $$ This follows from the fact that expectations are linear. Note here we assume $x_i$ is independent from $x_j$, unless of course $i = j$; this is how the expectation distributes over the product. If your inputs are correlated, this simple analysis fails. Note also that this is very close to what you would guess the "mean" value of $y$ is, but not quite. To second order, you have to account for the expectations of the squares of the inputs not matching the squares of the expectations. At least the maxim "plug in the best values of $x$ to get the best value of $y$" works perfectly to first order in $x$.

We can square this result, yielding \begin{align} \hat{y}^2 & = f_0^2 + 2f_0 \sum_{i=1}^n a_i \hat{x}_i + 4f_0 \sum_{i=1}^n \sum_{j=i+1}^n b_{ij} \hat{x}_i \hat{x}_j \\ & \quad \qquad + 2f_0 \sum_{i=1}^n b_{ii} E(x_i^2) + \sum_{i=1}^n \sum_{j=1}^n a_i a_j \hat{x}_i \hat{x}_j + \mathcal{O}(x^3). \end{align} Along similar lines, we have \begin{align} E(y^2) & = E\left(f_0^2 + 2f_0 \sum_{i=1}^n a_i x_i + 2f_0 \sum_{i=1}^n \sum_{j=i}^n b_{ij} x_i x_j + \sum_{i=1}^n \sum_{j=1}^n a_i a_j x_i x_j\right) + \mathcal{O}(x^3) \\ & = f_0^2 + 2f_0 \sum_{i=1}^n a_i \hat{x}_i + 4f_0 \sum_{i=1}^n \sum_{j=i+1}^n b_{ij} \hat{x}_i \hat{x}_j + 2f_0 \sum_{i=1}^n b_{ii} E(x_i^2) \\ & \quad \qquad + 2 \sum_{i=1}^n \sum_{j=i+1}^n a_i a_j \hat{x}_i \hat{x}_j + \sum_{i=1}^n a_i^2 E(x_i^2) + \mathcal{O}(x^3). \end{align}

Finally, we are in a position to compute the variance of $y$. This is simply $$ (\Delta y)^2 \equiv E(y^2) - \hat{y}^2 = \sum_{i=1}^n a_i^2 \left(E(x_i^2) - \hat{x}_i^2\right) + \mathcal{O}(x^3). $$ In fact, this result can be written entirely in terms of the first derivatives, $a_i$, and the variances of the inputs, $(\Delta x_i)^2$: $$ (\Delta y)^2 \approx \sum_{i=1}^n \left(\frac{\partial f}{\partial x_i}\right)^2 (\Delta x_i)^2. $$ The standard deviation $\Delta y$, which is what you expect reported when you see $y = \hat{y} \pm \Delta y$, is simply the square root of this: $$ \Delta y \approx \left(\sum_{i=1}^n \left(\frac{\partial f}{\partial x_i}\right)^2 (\Delta x_i)^2\right)^{1/2}. $$ It is this formula that people have in mind when they say things like "we added the uncertainties in quadrature," since this result is very broadly applicable - it applies whenever we have a continuously differentiable function of independent inputs with known means and standard deviations.

The idea is just that if the uncertainties are small enough you can approximate the function by its Taylor series

$$ f(x_i + \delta_i) \approx f(x_i) + \sum_j \frac{\partial f(x_i)}{\partial x_j} \delta_j + \sum_{j,k} \frac{1}{2} \frac{\partial^2 f(x_i)}{\partial x_j \partial x_k} \delta_j \delta_k + \cdots. $$

If you neglect the second order terms the function is a linear function and the uncertainties combine as they would for a linear function. If the uncertainties are so large that the higher order terms matter then you need to do something more complicated to propagate the uncertainty (e.g. Monte Carlo).