Why isn't my model learning satisfactorily? The problem to solve is non-linear regression of a non-linear function. My actual problem is to model the function "find the max over many quadratic forms": max(w.H.T * Q * w), but to get started and to learn more about neural networks, I created a toy example for a non-linear regression task, using Pytorch. 
The problem is that the network never learns the function in a satisfactory way, even though my model is quite large with multiple layers (see below). Or is it not large enough or too large? How can the network be improved or maybe even simplified to get a much smaller training error?

I experimented with different network architectures, but the result is never satisfactory. Usually, the error is quite small within the input interval around 0, but the network is not able to get good weights for the regions at the boundary of the interval (see plots below). The loss does not improve after a certain number of epochs.  I could generate even more training data, but I have not yet understood completely, how the training can be improved (tuning parameters such as batch size, amount of data, number of layers, normalizing input (output?) data,  number of neurons, epochs, etc.)

My neural network has 8 layers with the following number of neurons: 1, 80, 70, 60, 40, 40, 20, 1.

For the moment, I do not care too much about overfitting, my goal is to understand, why a certain network architecture/certain hyperparameters need to be chosen. Of course, avoiding overfitting at the same time would be a bonus.

I am especially interested in using neural networks for regression tasks or as function approximators. In principle, my problem should be able to be approximated to arbitrary accuracy by a single layer neural network, according to the universal approximation theorem, isn’t this correct?
Neural networks learn badly with large input ranges. Scale your inputs to a smaller range e.g. -2 to 2, and convert to/from this range to represent your function interval consistently.