I'm not an expert in neural nets but I think the following points might be helpful to you. There are also some nice posts, e.g. this one on hidden units, that you can search for on this site about what neural nets do that you might find useful.
1 Large errors: why didn't your example work at all
why errors are so big and why all predicted values are almost constant?
This is because the neural network was unable to compute the multiplication function you gave it and outputting a constant number in the middle of the range of y, regardless of x, was the best way to minimize errors during training. (Notice how 58749 is pretty close to the mean of multiplying two numbers between 1 and 500 together.)
It's very hard to see how a neural network could compute a multiplication function in a sensible way. Think about how each node in the network combines previously computed results: you take a weighted sum of the outputs from previous nodes (and then apply a sigmoidal function to it, see, e.g. an Introduction to Neural Networks, to scrunch the output inbetween $-1$ and $1$). How are you going to get a weighted sum to give you multiplication of two inputs? (I suppose, however, that it might be possible to take a large number of hidden layers to get multiplication working in a very contrived way.)
2 Local minima: why a theoretically reasonable example might not work
However, even trying to do addition you run into problems in your example: the network doesn't train successfully. I believe that this is because of a second problem: getting local minima during the training. In fact, for addition, using two layers of 5 hidden units is much too complicated to compute addition. A network with no hidden units trains perfectly well:
x <- cbind(runif(50, min=1, max=500), runif(50, min=1, max=500))
y <- x[, 1] + x[, 2]
train <- data.frame(x, y)
n <- names(train)
f <- as.formula(paste('y ~', paste(n[!n %in% 'y'], collapse = ' + ')))
net <- neuralnet(f, train, hidden = 0, threshold=0.01)
print(net) # Error 0.00000001893602844
Of course, you could transform your original problem into an addition problem by taking logs, but I don't think this is what you want, so onwards...
3 Number of training examples compared to number of parameters to estimate
So what would be a reasonable way to test your neural net with two layers of 5 hidden units as you originally had? Neural nets are often used for classification, so deciding whether $\mathbf{x}\cdot\mathbf{k} > c$ seemed a reasonable choice of problem. I used $\mathbf{k} = (1, 2, 3, 4, 5)$ and $c = 3750$. Notice that there are several parameters to be learnt.
In the code below I take a very similar approach to yours except that I train two neural nets, one with 50 examples from the training set, and one with 500.
library(neuralnet)
set.seed(1) # make results reproducible
N=500
x <- cbind(runif(N, min=1, max=500), runif(N, min=1, max=500), runif(N, min=1, max=500), runif(N, min=1, max=500), runif(N, min=1, max=500))
y <- ifelse(x[,1] + 2*x[,1] + 3*x[,1] + 4*x[,1] + 5*x[,1] > 3750, 1, 0)
trainSMALL <- data.frame(x[1:(N/10),], y=y[1:(N/10)])
trainALL <- data.frame(x, y)
n <- names(trainSMALL)
f <- as.formula(paste('y ~', paste(n[!n %in% 'y'], collapse = ' + ')))
netSMALL <- neuralnet(f, trainSMALL, hidden = c(5,5), threshold = 0.01)
netALL <- neuralnet(f, trainALL, hidden = c(5,5), threshold = 0.01)
print(netSMALL) # error 4.117671763
print(netALL) # error 0.009598461875
# get a sense of accuracy w.r.t small training set (in-sample)
cbind(y, compute(netSMALL,x)$net.result)[1:10,]
y
[1,] 1 0.587903899825
[2,] 0 0.001158500142
[3,] 1 0.587903899825
[4,] 0 0.001158500281
[5,] 0 -0.003770868805
[6,] 0 0.587903899825
[7,] 1 0.587903899825
[8,] 0 0.001158500142
[9,] 0 0.587903899825
[10,] 1 0.587903899825
# get a sense of accuracy w.r.t full training set (in-sample)
cbind(y, compute(netALL,x)$net.result)[1:10,]
y
[1,] 1 1.0003618092051
[2,] 0 -0.0025677656844
[3,] 1 0.9999590121059
[4,] 0 -0.0003835722682
[5,] 0 -0.0003835722682
[6,] 0 -0.0003835722199
[7,] 1 1.0003618092051
[8,] 0 -0.0025677656844
[9,] 0 -0.0003835722682
[10,] 1 1.0003618092051
It is apparent that the netALL does a lot better! Why is this? Take a look at what you get with a plot(netALL) command:
I make it 66 parameters that are estimated during training (5 inputs and 1 bias input to each of 11 nodes). You can't reliably estimate 66 parameters with 50 training examples. I suspect in this case you might be able to cut down on the number of parameters to estimate by cutting down on the number of units. And you can see from constructing a neural network to do addition that a simpler neural network may be less likely to run into problems during training.
But as a general rule in any machine learning (including linear regression) you want to have a lot more training examples than parameters to estimate.
