You don't need to use a double for loop here, as it goes through all possible combinations of elements of a and e_ (a(1) & e(1), a(1) & e(2), a(1) & e(3), and so on, total 16 combinations), which is not what you want.
A double for loop is not complicated but it is not required here. A single for loop is what you need.
Additionally, a simpler approach is to vectorize the code -
a = [6378137.000,6378137.000,6377397.155,6378136.000];
e_ = [0.0818191910428151,0.0818191908426216,0.0818191908426216,0.0818191065283643];
%Random value
fi = pi/3;
%Change the display format
format long g
M = (a.*(1-(e_).^2)) ./ (sqrt(1-(e_).^2) .* (sin(fi))^2).^3
Note the use of element-wise operators .*, ./ and .^
