We are trying to solve numerically a kink solution in matlab using the sine-gordon equation, therefore followed the steps according to the paper (1) listed below. The sine gordon equation is a non-linear pde and the numerical solution was clearly explained in the paper (1). However, the problem is that if we use the initial conditions f(x) and g(x) (page 6 of the paper), the solution will have a disturbance and is not as clean as figure 5.4a (p. 6) from the paper it self. To solve this we ‘doubled’ the size of the soliton by replacing ‘x’ in the initial conditions for ‘x/2’. Therefore my question: What is the explanition that dividing x by 2 instead of taking just x justifies? Perhaps we made a mistake in the numerical model?

Thanks for your time.

**Source**

- The paper used: https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ws2016-2017/num_methods_i/sinegordon.pdf
- There are a few references from source 1 to chapter 4: https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ws2016-2017/num_methods_i/wave.pdf

**Matlab code**

The changes of the inital conditions are cleary commented in the matlab code.

clear all; clc; close all;
%% Define variables
L = 40; %Choosing a bigger L can make it look nicer, also change it in the function at the bottom of the code
delta_x = 0.05;
delta_t = 0.01;
c = 0.2;
t_end = 90;
T = t_end/delta_t;
M = L/delta_x;
N = M+1;
x = linspace(-L,L,M+1);
t = linspace(0,t_end,T+1);
%end
%% Calculate first time row
u_0 = transpose(f(x));
gamma_1 = transpose(g(x));
beta_1 = sin(u_0);
alpha = delta_t/delta_x;
a = 1-alpha^2;
b = alpha^2/2;
c = b;
A = diag(a*ones(1,N)) + diag(b*ones(N-1,1),1) + diag(c*ones(1,N-1),-1);
A(1,2) = A(1,2) + b;
A(N,M) = A(N,M) + b;
B = 2*A;
%% Calculate second time row
u_1 = delta_t*gamma_1 + A*u_0 - (delta_t^2/2)*beta_1;
u = zeros(N,T);
u(:,1) = u_0;
u(:,2) = u_1;
%% Calculate the full matrix
beta = zeros(N,T);
beta(:,1) = beta_1;
for j = 2:(T)
beta(:,j) = sin(u(:,j));
u(:,j+1) = -u(:,j-1) + B*u(:,j) - (delta_t^2)*beta(:,j);
end
%figure
mesh(x(1:10:end),t(1:10:end),transpose(u(1:10:end,1:10:end)))
xlabel('x')
ylabel('t')
zlabel('u')
for n = 1:T-1
for l = 1:N-1
e(l,n) = (1/2*((u(l,n+1)-u(l,n))/delta_t)^2+1/2*((u(l+1,n+1)-u(l,n+1))/delta_x)*((u(l+1,n)-u(l,n))/delta_x)+(1-cos(u(l,n+1)+1-cos(u(l,n))))/2);
end
E(n) = delta_x*sum(e(:,n));
end
%% Define the functions
% Single soliton
function f = f(x)
c = 0.2; %This can be varied to create a nicer picture, also change it in the 'g'function
L = 20; %This can be varied to create a nicer picture, also change it in the 'g'function
f = 4*atan(exp((x/2)/sqrt(1-c^2))); % the '+L/2' can also be changed to have a different starting point for the soliton, don't forget to also change it in the 'g'function
end
function g = g(x)
c = 0.2;
L = 20;
g = -2*c/sqrt(1-c^2)*sech((x/2)/sqrt(1-c^2));
end