Using the GA R package to optimize the weights of a MLP neural network
I finally got the results I was looking for.

Apparently a single hidden layer neural network cannot learn the cubic function. Instead of just evolving the weights, I included the topology of the neural network into the chromosome. My neural network can have maximum 5 hidden layers with maximum 5 neurons in each layer and the chromosome contains information for the structure of the network and the selection of the weights. 

When I ran the GA with population of 200 and for 100 iterations I got a full topology of:

1(input) - 4-2-3-2-2 - 1(output), a mean squared error of 0.0051 and the following plot

I would be glad if anyone has any comments or questions.