Neural network cannot approximate complex nonlinear function

I have encountered a problem in training neural network. Even if I tried to increase hidden depth and width of hidden layers, the result turned out that it just cannot learn to approximate a nonlienear function. I enclose my code below.

n = 1;                           %input number
m = 1;                           %output number
hidden_depth = 3;                %hidden depth
hidden_width = 4;                %hidden width

Y_train_lags = linspace(-1, 1);
N_train = numel(Y_train_lags);
Y_train_diff = 16 * Y_train_lags.^5 - 20 * Y_train_lags.^3 + 5 * Y_train_lags;
Y_test_lags = linspace(-1, 1, 1000);
N_test = numel(Y_test_lags);
Y_test_diff = 16 * Y_test_lags.^5 - 20 * Y_test_lags.^3 + 5 * Y_test_lags;       %data used for training and test

layers = imageInputLayer([1 1 n]);
for i = 1:hidden_depth
    layers = [layers, fullyConnectedLayer(hidden_width), leakyReluLayer];
end
    layers = [layers, fullyConnectedLayer(m), regressionLayer];       %create layers with required depth and width

options = trainingOptions('adam', ...
    'MiniBatchSize', round(N_train/10), ...
    'LearnRateSchedule', 'piecewise', ...
    'LearnRateDropFactor', 0.9, ...
    'LearnRateDropPeriod', 20, ...
    'Shuffle', 'every-epoch', ...
    'Plots','training-progress', ...
    'MaxEpochs', 100, ...
    'InitialLearnRate', 1e-2);

net = trainNetwork(permute(Y_train_lags, [3 4 1 2]), Y_train_diff', layers, options);
Y_test_diff_hat = predict(net, permute(Y_test_lags, [3 4 1 2]))';

figure, plot(Y_train_lags, Y_train_diff, Y_test_lags, Y_test_diff_hat, '.');         %compare learning result and initial one

General result of running is as below.

Whatever how big depth and width I choose, the result of learning is still not good enough. I think the reason is the condition that all weights and bias equals to zeros is a bad local optimize for gradient descent. But I do not know how to let the program know what it converges to is a bad local extremum. And also, make the program to find a better extremum or even global optimization. Is there any existing function can do that? And is there any way I can change the initial value of weights and bias? It says that they are read-only.