Probably you need to do one or more of:
Decrease learning rate. Diverging loss is often a symptom of learning rate too high.
Increase number of hidden neurons. The output function will be the combination of many "patches" each created by a neuron that has learnt a different bias. The shape and quality of each patch is determined by the activation functions, but almost any nonlinear activation function used in an NN library should work to make a universal function approximator.
Normalise inputs. If you are training with high values for $x_1$ or $x_2$, this could make it harder to train unless you normalise your training data.
For your purposes, it might be an idea to skip the need for normalisation by training with $x_1$ and $x_2$ in range -2.0 to 2.0. It doesn't change your goal much to do this, and removes one potential problem.
You should note that a network trained against an open-ended function like this (where there are no logical bounds on the input params) will not learn to extrapolate. It never "learns" the function itself, but an approximation of it close to the supplied examples. When you supply a $(x_1, x_2)$ input far from the training examples, the output is likely to completely mismatch your original function.