Source Localization with the BLUE Estimator
One problem with your code is that the range calculations should be

R(i) = sqrt((xs-x(i))^2$\color{red}{\bf +}$(ys-y(i))^2);

Rn(i) = sqrt((xn-x(i))^2$\color{red}{\bf +}$(yn-y(i))^2);

But that doesn't appear to help. 

The thing that seems to be the issue is that your noise variance is way too high.  If I set it to var = 0.00000001 then I get good estimates.

The reason is because you're using this variance (or a standard deviation of 0.0001) to add noise to your H matrix. Looking at the values there:



you can see you're adding something with a variance of 1 (standard deviation of 1), to something of the order of $10^{-5}$.  That means your SNR (signal-to-noise ratio) is about -100dB.  I know of no algorithm that will help you at SNRs that low.



You left out that Kay defines
$$
R_i(x_s,y_s) = \sqrt{(x_s - x_i)^2 + (y_s - y_i)^2}
$$
and
$$
R_{n_i}(x_n,y_n) = \sqrt{(x_n - x_i)^2 + (y_n - y_i)^2}
$$
where $x_s = x_n + \delta x_s$ and $y_s = y_n + \delta y_s$.

Then the aim is to write $R_i$ in terms of $R_{n_i}$.

Wikipedia has an example of Taylor series in two variables, and the first order approximation is just:

$$
f(x,y) \approx f(a,b) + (x-a) \frac{df}{dx} (a,b) + (y-b) \frac{df}{dy} (a,b) + \ldots
$$

Then all you need to do is apply Wikipedia's expression to Kay's to get

$$
R_i \approx R_{n_i} + \frac{x_n-x_i}{R_{n_i}} \delta x_s + \frac{y_n - y_i}{R_{n_i}} \delta y_s. 
$$