Backlund Transform for multi-solitons of sine-Gordon model

p1[y_]:=4 ArcTan[Exp[y]];
ys[x_,t_,a_]:=((a+1/a) x+(a-1/a) t)/2; (* gamma (x-v t) v= (1/a-a)/(1/a+a) *)
p2[y1_,y2_,a1_,a2_]:=4 ArcTan[(a1+a2)/(a1-a2) Tan[(p1[y1]-p1[y2])/4]];
p3[y1_,y2_,y3_,a1_,a2_,a3_]:=p1[y1]+4 ArcTan[(a2+a3)/(a2-a3) Tan[(p2[y1,y2,a1,a2]-p2[y1,y3,a1,a3])/4]];
p4[y1_,y2_,y3_,y4_,a1_,a2_,a3_,a4_]:=p2[y1,y2,a1,a2]+4 ArcTan[(a3+a4)/(a3-a4) Tan[(p3[y1,y2,y3,a1,a2,a3]-p3[y1,y2,y4,a1,a2,a4])/4]];

Plot[Log[v]/(v/Sqrt[1-v^2]),{v,.1,.9}]

a1=1.2;a2=1/a1;a3=.01; t1=217;
Animate[Plot[p1[ys[x+40,t,a2]],{x,-60,60},PlotRange->{-8,15}],{t,0,t1}]
Animate[Plot[p2[ys[x+40,t,a2],ys[x-40,t,a1],a1,a2],{x,-60,60},PlotRange->{-8,15}],{t,0,t1}]

Energy Plot

one-soliton

Simplify[D[p1[ys[x,t,a]],x]^2+D[p1[ys[x,t,a]],t]^2]

Energy1[x_,t_,a_]:=1/2 (((8 (1+(a^4)) (E^(((1+(a^2)) (t+x))/a)))/((a^2) ((((E^(t/a))+(E^((x/a)+a (t+x)))))^2))))+(1-Cos[p1[ys[x,t,a]]]);
a11=.5;
v[a_]:=(1/a-a)/(1/a+a);
v[a11]
Animate[Plot[Energy1[x+10,t,a11],{x,-10,10},PlotRange->{0,8}],{t,0,40}]

Plot[{Energy1[x,0,0.5],Energy1[x,0,0.7],Energy1[x,0,0.9]},{x,-5,5},PlotRange->{0,7}]
v[0.5]
v[0.7]
v[0.9]

two-soliton

p2x[x_,t_,a1_,a2_,x1_,x2_]:=D[p2[ys[x-x1,t,a1],ys[x-x2,t,a2],a1,a2],x];
p2t[x_,t_,a1_,a2_,x1_,x2_]:=D[p2[ys[x-x1,t,a1],ys[x-x2,t,a2],a1,a2],t];
kineticE2[x_,t_,a_,b_,c_,d_]:=1/2((((a+b))^2) ((((((2 (-1+(a^2)) (E^((c-(a^2) c+(1+(a^2)) (t+x))/(2 a))))/(a ((E^((c+t)/a))+(E^((x/a)+a (-c+t+x))))))-((2 (-1+(b^2)) (E^((d-(b^2) d+(1+(b^2)) (t+x))/(2 b))))/(b ((E^((d+t)/b))+(E^((x/b)+b (-d+t+x))))))))^2)+(((((2 (1+(a^2)) (E^((c-(a^2) c+(1+(a^2)) (t+x))/(2 a))))/(a ((E^((c+t)/a))+(E^((x/a)+a (-c+t+x))))))-((2 (1+(b^2)) (E^((d-(b^2) d+(1+(b^2)) (t+x))/(2 b))))/(b ((E^((d+t)/b))+(E^((x/b)+b (-d+t+x))))))))^2)) ((Sec[ArcTan[(E^((1/2) (-((c+t-x)/a)+a (-c+t+x))))]-ArcTan[(E^((1/2) (-((d+t-x)/b)+b (-d+t+x))))]])^4))/((((a-b))^2) (((1+(1/(((a-b))^2))(((a+b))^2) ((Tan[ArcTan[(E^((1/2) (-((c+t-x)/a)+a (-c+t+x))))]-ArcTan[(E^((1/2) (-((d+t-x)/b)+b (-d+t+x))))]])^2)))^2));
Energy2[x_,t_,a1_,a2_,x1_,x2_]:=kineticE2[x,t,a1,a2,x1,x2]+(1-Cos[p2[ys[x-x1,t,a1],ys[x-x2,t,a2],a1,a2]]);

Needs["NumericalCalculus`"]

NEnergy2[x_,t_,a_,b_,c_,d_]:=1/2 (ND[p2[ys[x-c,tc,a],ys[x-d,tc,b],a,b],tc,t])^2+1/2 (ND[p2[ys[xc-c,t,a],ys[xc-d,t,b],a,b],xc,x])^2+(1-Cos[p2[ys[x-c,t,a],ys[x-d,t,b],a,b]]);
NEnergy3[x_,t_,a_,b_,c_,d_,e_,f_]:=1/2 (ND[p3[ys[x-d,tc,a],ys[x-e,tc,b],ys[x-f,tc,c],a,b,c],tc,t])^2+1/2 (ND[p3[ys[xc-d,t,a],ys[xc-e,t,b],ys[xc-f,t,c],a,b,c],xc,x])^2+(1-Cos[p3[ys[x-d,t,a],ys[x-e,t,b],ys[x-f,t,c],a,b,c]]);

a2=.6;a1=1;a3=.5;x1=30.;x2=-30.;x3=-1.;Tmax=200;
v[z_]:=(1/z-z)/(1/z+z);
x2f=v[a2] Tmax+x2
Animate[Plot[Energy2[x,t,a1,a2,x1,x2],{x,-50,50},PlotRange->{0,10}],{t,0,Tmax}]

Plot[NEnergy2[x,0,a1,a2,x1,x2],{x,-12,50},PlotRange->{0,10}]

Animate[Plot[NEnergy2[x,t,a1,a2,x1,x2],{x,-12,50},PlotRange->{0,10}],{t,0,Tmax}]

three-soliton

Animate[Plot[NEnergy3[x,t,a1,a2,a3,x1,x2,x3],{x,-12,50},PlotRange->{0,20}],{t,0,Tmax}]

b2=.7;b1=1;y1=0.;y2=-10.;
y2f=v[b2] Tmax+y2
Animate[Plot[Energy2[x,t,b1,b2,y1,y2],{x,-12,50},PlotRange->{0,6}],{t,0,Tmax}]