%根據k1和k2的值計算雙孤子解v的表達式
%
%

clear all;

%k1=10;
%k2=-10;

v='-1/2*(4*exp(2*z1+3*z2+A)*k1^2+3*exp(3*z1+2*z2+A)*k1^2+2*exp(3*z1+z2+A)*k2^2+k1^2*exp(3*z1)+k2^2*exp(3*z2)+4*k2^2*exp(2*z1+2*z2)-4*k1*exp(2*z1+z2)*k2+2*k2^2*exp(z1+3*z2)-8*k1*exp(2*z1+2*z2)*k2-4*k1*exp(z1+2*z2)*k2-2*exp(3*z1+3*z2+2*A)*k1*k2+2*k1^2*exp(z1+3*z2)-4*exp(3*z1+2*z2+A)*k1*k2-4*k1*exp(z1+3*z2)*k2+4*k1^2*exp(2*z1+2*z2)+4*exp(3*z1+2*z2+A)*k2^2-4*exp(2*z1+3*z2+A)*k1*k2+3*exp(2*z1+3*z2+A)*k2^2+4*k2^2*exp(z1+2*z2)+2*exp(z1+3*z2+A)*k2^2+3*k2^2*exp(2*z1+z2)+exp(3*z1+3*z2+2*A)*k1^2+exp(4*z1+z2+A)*k2^2+3*k1^2*exp(z1+2*z2)+exp(3*z1+3*z2+2*A)*k2^2+2*exp(3*z1+z2+A)*k1^2+exp(z1+4*z2+A)*k1^2+2*k2^2*exp(3*z1+z2)+exp(3*z1+2*z2+2*A)*k1^2+4*k1^2*exp(2*z1+z2)+2*k1^2*exp(3*z1+z2)+4*exp(2*z1+2*z2+A)*k1^2+4*exp(2*z1+2*z2+A)*k2^2+exp(2*z1+z2+A)*k2^2+exp(z1+2*z2+A)*k1^2+k1^2*exp(z1+z2)+k2^2*exp(z1+z2)+2*exp(z1+3*z2+A)*k1^2-2*k1*exp(z1+z2)*k2-4*k1*exp(3*z1+z2)*k2+exp(2*z1+3*z2+2*A)*k2^2)/(exp(z1)+exp(z2)+exp(z1+z2+A))^2/(1+exp(z1)+exp(z2))^2';

v=subs(v,'k1','k');
v=subs(v,'k2','k');

v=subs(v,'2*A','');  %注意先替換'2*A'
v=subs(v,'A','');

v=subs(v,'4*z1','q4');  %注意替換的順序
v=subs(v,'3*z1','q3');
v=subs(v,'2*z1','q2');
v=subs(v,'z1','q1');    %注意這里沒有'1*'

v=subs(v,'4*z2','r4');  %注意替換的順序
v=subs(v,'3*z2','r3');
v=subs(v,'2*z2','r2');
v=subs(v,'z2','r1');    %注意這里沒有'1*'

v=simplify(v)

%后面的化簡需要人工介入
%在Matlab中顯示v，拷貝v的公式到記事本中，后面的操作在記事本中完成
%去掉分母，因為分母單獨計算
%將指數里面的'+'換為'.*'，即點乘運算，不屬于指數里面的'+'不要替換
%將‘exp'換為空，即去掉指數，因為在上一步中已經指數運算化成了點乘運算
%用s點乘v，即v=s.*v;
%拷貝公式到繪圖程序中的適當位置用于計算v





%當k1=k2=k時的表達式
%v ='-1/2*k^2*(3*exp(q3+r2+qa)+3*exp(q1+r2)+exp(q3+r2+q2a)+exp(q1+r2+qa)+8*exp(q2+r2+qa)+3*exp(q2+r3+qa)+exp(q2+r3+q2a)+4*exp(q1+r3+qa)+exp(q4+r1+qa)+exp(q2+r1+qa)+4*exp(q3+r1+qa)+3*exp(q2+r1)+exp(q3)+exp(r3)+exp(q1+r4+qa))/(exp(q1)+exp(r1)+exp(q1+r1+qa))^2/(1+exp(q1)+exp(r1))^2'


%當k1=k,k2=-k時的表達式
%v
%='-1/2*k^2*(11*exp(q2+r3+qa)+8*exp(q3+r1)+exp(q4+r1+qa)+exp(q2+r1+qa)+4*exp(q1+r1)+4*exp(q3+r1+qa)+11*exp(q2+r1)+exp(q2+r3+q2a)+4*exp(q1+r3+qa)+8*exp(q1+r3)+4*exp(q3+r3+q2a)+exp(q1+r4+qa)+exp(r3)+8*exp(q2+r2+qa)+exp(q3)+exp(q3+r2+q2a)+exp(q1+r2+qa)+11*exp(q3+r2+qa)+16*exp(q2+r2)+11*exp(q1+r2))/(exp(q1)+exp(r1)+exp(q1+r1+qa))^2/(1+exp(q1)+exp(r1))^2'