How to numerically estimate MLE estimators in python when gradients are very small far from the optimal solution?
There are several reasons why you are getting erroneous results.
First, you should consider using log likelihood instead of likelihood. There are numerical issues with multiplying many small numbers(imagine if you had millions of samples you had to multiply millions of small numbers for the lhd). Also taking gradients for optimization methods that require gradients is often easier when you are dealing with the log likelihood. In general, it is good to have an objective which is a sum rather than a product of variables when dealing with optimization problems.
Second, fmin is using Nelder-Mead simplex algorithm which has no convergence guarantees according to scipy documentation. This means the convergence is totally random and you should not expect to find parameters close to the originals. To get around this, I would suggest you to use a gradient based method like stochastic gradient descent or BFGS. Since you know the generative model (rvs are Gaussian distributed) you can write the likelihood and log likelihood as:

Where a,b,c and d are your model parameters 5,2,3 and 4 respectively.
Then take the gradient with respect to [a,b,c,d] and feed that into prime input of fmin_bfgs. Note that due to varying variance what could be solved by just linear regression is now a nastier problem.
Finally, you may also want to check Generalized least squares here and here, which talk about your problem and offer several available solutions.
Good luck!