function myTSP2
%城市數目
N=10;
%城市坐標及城市間距離
cityx=[0.4,0.2439,0.1707,0.2293,0.5171,0.8732,0.6878,0.8488,0.6683,0.6195];
cityy=[0.4439,0.1463,0.2293,0.761,0.9414,0.6536,0.5219,0.3609,0.2536,0.2634];
for i=1:1:10
    for j=1:1:10
        d(i,j)=sqrt((cityx(i)-cityx(j))^2+(cityy(i)-cityy(j))^2);
    end
end
%網絡參數
A=500;
B=500;
C=200;
D=500;
u0=0.02;
tao=1;
lamda=0.00001;
%多次實驗并統計結果
for count=1:1:10
    %換位陣及初始化
    V=0.5+rand(N,N)*0.1;
    U=atanh(2*V-1)*u0;
    %狀態更新
    E1=0;
    toend=0;
    while toend==0;
        %同步更新
        for ux=1:1:N
            for ui=1:1:N
                m1=0;
                m2=0;
                m3=0;
                m4=0;
                %求導公式第一項
                for j=1:1:N
                    if j~=ui
                        m1=m1+V(ux,j);
                    end                   
                end
                m1=-A*m1;
                %求導公式第二項
                for y=1:1:N
                    if y~=ux
                        m2=m2+V(y,ui);
                    end
                end
                m2=-B*m2;
                %求導公式第三項
                for x=1:1:N
                    for j=1:1:N
                        m3=m3+V(x,j);
                    end
                end
                m3=-C*(m3-N);
                %求導公式第四項
                for y=1:1:N
                    if y~=ux && ui<N && ui>1
                        m4=m4+d(ux,y)*(V(y,ui+1)+V(y,ui-1));
                    end
                    if y~=ux && ui==1
                        m4=m4+d(ux,y)*(V(y,ui+1)+V(y,N));
                    end
                    if y~=ux && ui==N
                        m4=m4+d(ux,y)*(V(y,ui-1)+V(y,1));
                    end
                end
                m4=-D*m4;
                Udao(ux,ui)=-U(ux,ui)+m1+m2+m3+m4;
            end
        end
        %導數及狀態更新                
        U=U+lamda*Udao;
        V=(1+tanh(U/u0))/2;     
        for ux=1:1:N
            for ui=1:1:N
                if V(ux,ui)<0.01
                    V(ux,ui)=0;
                end
                if V(ux,ui)>0.9
                    V(ux,ui)=1;
                end 
            end
        end
        V ;
        %求能量函數
        e1=0;
        e2=0;
        e3=0;
        e4=0;
        for x=1:1:N
            for i=1:1:N-1
                for j=i+1:1:N
                    e1=e1+V(x,i)*V(x,j);
                end
            end
        end
        e1=A*e1;
        for i=1:1:N
            for x=1:1:N-1
                for y=x+1:1:N
                    e2=e2+V(x,i)*V(y,i);
                end
            end
        end
        for i=1:1:N
            for i=1:1:N
                e3=e3+V(x,i)*V(x,i);
            end
        end
        e3=C*(e3-N)^2;
        for x=1:1:N
            for y=1:1:N
                for i=1:1:N
                    if i==1
                        e4=e4+d(x,y)*V(x,i)*(V(y,i+1)+V(y,N));
                    end
                    if i==N
                         e4=e4+d(x,y)*V(x,i)*(V(y,1)+V(y,i-1));
                    end
                    if i~=1 && i~=N
                        e4=e4+d(x,y)*V(x,i)*(V(y,i+1)+V(y,i-1));
                    end
                end
            end
        end
        e4=D*e4;
        E0=E1;
        E1=e1+e2+e3+e4;    
        if E1==E0%(E1-E0)^2<0.00000001
            toend=1;
        else
            toend=0;
        end         
    end
    V
end
V
for j=1:1:N
    for i=1:1:N
        if V(i,j)==1
            cityx_final(j)=cityx(i);
            cityy_final(j)=cityy(j);
        end
    end
end
cityx_final(N+1)=cityx(1);
cityy_final(N+1)=cityy(1);
cityx_final
cityy_final
plot(cityy_final,cityx_final);