Positive Index k=1010

function xy = solver(kd,bx)

p0=bx([1 3]);
npts = length(kd);
nones=ones(npts,1);
if any(kd(:)>0)
  dx=bx(2)-bx(1);
  dy=bx(4)-bx(3);
D=dx+dy;
d=D/30;
while max(2,ceil(dx/d))*max(2,ceil(dy/d))<200
	d=d/2;
end
x=0:d:dx;
if x(end)<dx
	x(end+1)=dx;
end
y=0:d:dy;
if y(end)<dy
	y(end+1)=dy;
end
[x,y]=meshgrid(x,y);
xy0=x(:).'+i*y(:).';
k = length(xy0);
  [dum,p]=sort(-max(kd));
  kd=kd(p,p);
  S=(kd>0); % Remember positive position
  sI=find(S); % Positive index
  strainMatrix = S; % Create strainMatrix;
  M=9;
  N0=min(M,npts);
  xy = zeros(npts,1);
  xy(1:N0) = solverl(S(1:N0,1:N0),kd(1:N0,1:N0),dx,dy);
  for N=(M+1):npts
    xy(N) = solver2(S(1:N,1:N),kd(1:N,1:N),N,xy0,...
             xy(1:N,:),k,d,dx,dy);
  end
  for N=1:N0
    xy(N) = solver2(S,kd,N,xy0,xy,k,d,dx,dy);
  end
  i0=0;
  for N=1:npts
    [XY1,XY2] = meshgrid(xy,xy.');
    strainMatrix(sI) = abs(abs(XY1(sI)-XY2(sI))-kd(sI));
    [m,i]=max(sum(strainMatrix));
    if i==i0, break; else, i0=i; end
    xy(i) = solver2(S,kd,i,xy0,xy,k,d,dx,dy);
  end
  [v,w]=sort(p);
  xy=xy(w,:);
  xy=[real(xy) imag(xy)]+p0(nones,:);
else
  xy = p0(nones,:);
end


function xy = solverl(S,kd,dx,dy)
npts = size(S,1);
[I,J] = find(triu(S));
y = kd(find(triu(S)));
k = 1010; % change k according to memory capacity.
x = rand(npts,k)*dx+i*rand(npts,k)*dy;
[tv,id] = min(sum(abs(abs(x(I,:)-x(J,:))-y(:,ones(1,k)))));
xy=x(:,id);


function xy = solver2(S,kd,N,xy,xy0,k,d,dx,dy)
I = find(S(:,N));
if isempty(I),
  xy = xy0(N);
  return;
end
z = kd(:,N);
kones=ones(1,k);
[XY1,XY2]=meshgrid(xy,xy0(I));
[tv,id]=min(sum(abs(abs(XY1-XY2)-z(I,kones))));
x1=real(xy(id));
y1=imag(xy(id));
x=max(0,x1-d):d/10:min(dx,x1+d);
y=max(0,y1-d):d/10:min(dy,y1+d);
[x,y]=meshgrid(x,y);
xy=(x+i*y);
xy=xy(:).';
k = length(xy);
kones=ones(1,k);
[XY1,XY2]=meshgrid(xy,xy0(I));
[tv,id]=min(sum(abs(abs(XY1-XY2)-z(I,kones))));
xy=xy(id);