@Thomas Dean Yes, plot and plot3d commands (and maybe a few other plotting commands) use, when possible, something called evalhf instead of evalf. The evalhf tends to be 20-30 times faster. The h stands for hardware.

Here's how evalhf can be applied to the current example:

restart:
st:= time():
y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m):
f:= unapply(abs(y), m):
n := 500: ## sample size
M := <seq(2.*idx/n,idx=1..n)>: ## m
Y := map[evalhf](f,M)+~Statistics:-Sample(Normal(0,3), n)^+: ## signal + noise
time()-st;

restart:
y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m):
f:= (x) -> abs(subs(m=x,y)):
n := 500: ## sample size
M := <seq(2.*idx/n,idx=1..n)>: ## m
Y := map[evalhf](f,M)+~Statistics:-Sample(Normal(0,3), n)^+: ## signal + noise

                             0.062
Error, unable to evaluate built-in function `subs` in evalh

See? The evalhf is very restricted in the commands that it can use. They're listed at ?evalhf,fnclist.

@samira moradi Use command Export, ExportMatrix, or ExcelTools:-Export. They all get the job done one way or another.

@Thomas Dean The notation A:-B is definitely superior, more robust, and clearer to read, and less ambiguous than A[B]. It's a shame that the latter is emphasized in the help pages. There is a rare situation where A[B] is required: When B isn't literally the name of an export of module A but is an expression which evaluates to the name of an export.

I took your code, commented out most lines, and wrote under those lines my suggested improvements:

restart:
with(Statistics):
y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m);
## fix complex values
#f := (x) -> abs(evalf(subs(m=x,y)));
f:= unapply(abs(y), m);

n := 500;

## mean 1, std dev 0.5
#X := RandomVariable(Normal(1, .5)):
M := Sample(Normal(1,.5), n):
#Y := [seq(f(M[idx]),idx=1..n)]:
Y:= f~(M):

#A:=[seq(idx,idx=1..n)]:
#plot(A,M,style=point,symbol= circle, symbolsize= 1,title="Values of M");
Statistics:-Histogram(M, title= "Values of M", frequencyscale= absolute);
#plot(A,Y,style=point,symbol= circle, symbolsize= 1,title="Values of Y");
Statistics:-Histogram(Y, title= "Values of Y", frequencyscale= absolute);
plot(Y,M,style=point,symbol= circle, symbolsize= 1,title="M vs Y");
#DataSummary(M); DataSummary(Y);
DataSummary~([M,Y])[];
#Data:=Matrix(n,2,[seq([Y[idx],M[idx]][],idx=1..n)]):
Data:= <Y,M>^+:

plot(
     [[y, m, m= 0..2], Data], labels= ['y', m],
     style= [line, point], color= [red, black],
     symbol= point, symbolsize= 1
);

@Markiyan Hirnyk Sorry, you're right, I used the word hyperplane incorrectly. I meant a nontrivial affine subspace.

@Thomas Dean To generate the Data, I extracted the point data from the first plot and added a Normal(0, .1) to each second coordinate. Like this:

y := int(1/(-0.4016e-1*m^(2/3)-0.211e-3*m^(5/3)), m):
P:= plot([y, m, m= 0..2]):
Data:= op([1,1], P):
Data[..,2]:= Data[..,2] + Statistics:-Sample(Normal(0,.1), op([1,1], Data))^+:

That's all there is to it.

@Ronan See ?use, ?object,operators, ?ProgrammingGuide,chapter09.

@Markiyan Hirnyk Agreed. In some ways Explore is better than an animation (flexibility). The animation is automatic and exportable.

@vv Iterator overwrites the same Array on each iteration rather than generating a new Array. So it seems likely that the indices would be returned in the same order on each iteration. Only a kernel programmer could say for sure.

@Markiyan Hirnyk Both techniques give solutions that minimize the 2-norm of the residuals. Any differences there are due solely to round-off error and can be resolved by manipulating Digits. But there's also the issue of minimizing the 2-norm of the results (the C[i] values) themselves. That's what optimize is doing in this case. This system is unusual in that it's both inconsistent and there's a hyperplane of residual-minimizing solutions.

@samira moradi Suppose that Data is an n x 2 Matrix containing the data. In this case, the y data should be in the first column, and the m data should be in the second column. Then simply modify the first plot command to

plot(
     [[y, m, m= 0..2], Data], labels= ['y', m],
     style= [line, point], color= [red, black],
     symbol= point, symbolsize= 1
);

@Markiyan Hirnyk My Answer is different and better and was written without having seen your Answer. It can be expressed in a single line, and it gives the solution of minimal 2-norm.

@Markiyan Hirnyk DirectSearch has returned some very large values, although correct. LeastSquares with the optimize option provides the minimal-2-norm solution. 

@vv I saw that the usage on Arrays (of any dimension) was documented and unfortunately deleted my Reply before you Replied.

@vv I think that your usage of entries is risky if you don't include option indexorder. It may be that the entries of a Vector always come in index order, but I can't find that documented. Since it's a kernel command, we can't look at its code.

Nonetheless, vote up.

@Joe Riel Aha! I've encountered this before. The problem is the plus sign!