From a lecture in machine learning, I know that a linear activation function can only produce a linear function, but I don't know if it can produce a connected linear function, like the one in the image below. This function consists of multiple concatenated lines.
A linear activation would not be able to separate the data like you have shown, no matter how many layers you throw into the network.
If we had multiple linearly activated layers, each feeding into each other, the neurons in the previous layer would calculate some weighted sum of the input and send it to the next layer as input, where the next layer's neurons would also calculate a weighted sum on that input, and, in turn, it fires based on another linear activation function.
No matter how many layers and neurons there are, if all are linear in nature, the final activation function of the last layer is also a linear function of the input of the first layer.
That means that any or all of these layers can be replaced by one layer. This completely loses the advantage of stacking layers because any multilayer network is equivalent to a single layer with linear activation, given that a combination of linear functions in a linear manner is still another linear function.
A good way to visualize this messing around with Tensorflow playground, which has a spiral data set similar to your data.
In contrast to the failure of linear activation functions for this data in the previous link, check out this much smaller network using a Tanh activation function which can separate the data within 100 or so iterations.
For further reference on activation functions, checkout this blog post.
Since I can't comment, there are a few caveats to previous answers.
For instance, if you knew beforehand what the expected boundary function for that variable was, then you could transform it first. For instance, if you knew one feature was expected to be sinusoidal, you could transform your data (theta) using $f(x) = a*sin(\theta)$ first then expect the variable to be linear. At that point it would fit a linear function and solve for the amplitude instead. Of course, in large problems, this may be impractical, but still possible.
Likewise with quadratic data. I could add 1 new feature for every continuous variable that looks like $x^2$ and capture any quadratic dependence.
Another solution in your shown example would be if you knew the data was piece-wise, you could train 3 small, easily trained, NNs for instance with different boundaries. During each new training, you set new boundaries for your functions and see which produces the best results. In your case, this would be something close to $[-\infty,-0.5], [-0.5,0.5]$, and $[0.5,\infty]$. You could argue that this is a single piece-wise neural net just like any traditional piece-wise function. While this last approach seems very simplistic/impractical, this is actually done in practice. For instance in autonomous driving, a separate NN could be trained for object detection (human, lights, signs etc.), another for steering wheel positioning for different speeds,conditions, and turning radius, another for how to accelerate, etc. and then the system runs them all and combines their outputs into a full driving car. Each output is then its boundary of the puzzle on what region of "driving space" it handles.
So while I would say the simple answer is no, there are still ways around it if you have some prior knowledge.
As a response for the Universal Approximation Theorem, see this wiki, which says
One of the first versions of the theorem was proved by George Cybenko in 1989 for sigmoid activation functions. It was later shown that the class of deep neural networks is a universal approximator if and only if the activation function is not polynomial."