T=ones(17);
 Q=100;
x=[10 10 10 10 20 20 20 20 30 30 30 30 40 40 40 40];
y=[10 20 30 40 40 30 20 10 10 20 30 40 40 30 20 10];
axis([0 50 0 50 ])
 hold on
 plot(x,y,'*')
 hold on
for x1=1:16
      for y1=1:16
         d(x1,y1)=sqrt((x(x1)-x(y1))^2+(y(x1)-y(y1))^2);
      end
 end 
for i=1:20
  for b=1:20       
    v=1;
    p(:,17)=0;
    p(17,:)=0;
    m(1)=1;
    h=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17];
    q=0;
  for c=1:15              
    for k=1:16              
        if h(k)==m(v);
           h(k)=17;
           v=v+1;
           break;
        end
    end   
        for i1=1:16       
           for j1=1:16
               if   d(i1,j1)==0
                    z(i1,j1)=0;    
               else z(i1,j1)=1/d(i1,j1);
               end         
           end
        end
        z(:,17)=0;
        for j=1:16         
            p(m(v-1),h(j))=(T(m(v-1),h(j))^1)*(z(m(v-1),h(j))^5);
        end 
            p(m(v-1),17)=0;
        u=0;
        h1=h;
        h2=0;
        for k=1:16
            h0=m(v-1)-h1(k)
            if  abs(h0)>3;
                h1(k)=17;
            end
        end  
                for j=1:16
                     if h1(j)==17
                        h2=h2+1; 
                     end
                end
       if h2==16
          break;
       end       
        for k=1:16
            u=u+p(m(v-1),h1(k));
        end 

            
   t=-10; 
    for k=1:16               
        if p(m(v-1),h1(k))/u>t    
           t=p(m(v-1),h1(k))/u;
           m(v)=h(k);
        end
    end
if v<16
   if p(m(v-1),m(v))/u<rand(1)
      g=p(m(v-1),m(v));
      p(m(v-1),m(v))=0;
      t=-10;
      for k=1:16               
        if p(m(v-1),h1(k))/u>t    
           t=p(m(v-1),h1(k))/u;
           m(v)=h1(k);
        end
      end
      p(m(v-1),m(v))=g;
   end   
 end

end
 
w=0;
for k=1:16
    if m(k)==16
       w=w+1; 
    end 
   if h(k)~=17&w~=1
      L((i-1)*16+b)=100000;
   else   L((i-1)*16+b)=0;
     for s=1:15          
        L((i-1)*16+b)=L((i-1)*16+b)+sqrt((x(m(s))-x(m(s+1)))^2+(y(m(s))-y(m(s+1)))^2);
     end 
      
   end
end 
   for k=1:((i-1)*16+b)
    if w~=1 
       break;
    end
    if min(L)==L((i-1)*16+b)
       e=m;
       for s=1:15 
           T(m(s),m(s+1))=(1-0.99)*T(m(s),m(s+1))+0.01*Q/L((i-1)*16+b);
       end
      
     end
   end    
 end
end
  for r=1:15
       line([x(e(r)),x(e(r+1))],[y(e(r)),y(e(r+1))]);
  end