Why is non-linearity desirable in a neural network?
Consider what happens if you intend to train a linear classifier on replicating something trivial as the XOR function. If you program/train the classifier (of arbitrary size) such that it outputs XOR condition is met whenever feature a or feature b are present, then the linear classifier will also (incorrectly) output XOR condition is met whenever both features together are present. That is because linear classifiers simply sum up contributions of all features and work with the total weighted inputs they receive. For our example, that means that when the weighted contribution of either feature is sufficient already to trigger the classifier to output XOR condition is met, then obviously also the summed contributions of both features are sufficient to trigger the same response.
To get a classifier that is capable of outputting XOR condition is met if and only if the summed contributions of all input features are above a lower threshold and below an upper threshold, commonly non-linearities are introduced. You could of course also try to employ a quadratic function to solve the two-feature problem, but as soon as the number of variables/features increases again, you run into the same problem again, only in higher dimensions. Therefore, the most general approach to solving this problem of learning non-linear functions, like the XOR, is by setting up large models with enough capacity to learn a given task and equipping them with non-linearities. That simplifies training since it allows for using stochastic gradient descent for training the system/classifier, preventing one from having to solve Higher-Degree polynomial equations analytically (which can get computationally quite expensive quite quickly) to solve some task.
In case you are interested, here's one paper analyzing and dealing with the XOR problem (as one concrete instance of a problem where purely linear models fail to solve some task).
EDIT:
You can consider a layer in a network as a function $y = f(x)$, where $x$ is the input to some layer $f$ and $y$ is the output of the layer. As you propagate $x$, being the network's input, through the network, you get something like $y = p(t(h(g(f(x)))))$, where $f$ is the input layer and $p$ constitutes the output layer, i.e. a set of weights, by which the input to that respective layer gets multiplied.
If $h$, for example, is some non-linear activation function, like ReLU or sigmoid, then $y$, being the network's output, is a non-linear function of input $x$.