x = [535 749 963 1177 1391 1605 1819 2033 2247 2461 2675 2889]; %Distance from edge (m)
y = [-15.0 -24.0 -31.2 -36.8 -40.8 -42.7 -42.4 -40.9 -37.3 -31.5 -21.8 -12.8].*10.^-5; %Anomaly (mgal)
M = 1:34.2:3420; %M regular elements
dx = 34.2; %M is dx wide
%Plot data
figure(1);
plot(x,y,'*');
xlabel('Distance from west edge (m)');
ylabel('Bouguer anomaly (mgal)');
%The inverse problem
%(1) Perform discretization
G = 6.67*10^(-11); %Gravity constant in [m3/ kg*s^2]
del_rho = -1700; %Density contrast
d = 0.1;
h=zeros(1,1200);
for j = 1:12 %We have 12 measurements
for i = 1:100 %100 M regular elements
h((j-1)*100+i)=sqrt(exp(y(j)/(G*del_rho*dx))*(M(i)-x(j))^2 + x(j)^2 - M(i));%you can directly solve for h but there are 1200 different values
end
end
semilogy(1:1200,h);%there are 100 values of h for each of the 12 measurements, you could reshape(h,12,100) to get rows for each measurement.
I am not sure this is what you want, but you can easily solve for h in your equation.
