% SH-Modell mit Crank Nicholson periodische Randbedingungen
nwave=20
xmax=nwave*2*pi
alpha=.6
R=.5
fig1=figure(1);

nx=200
dx=xmax/(nx+1)
tmax=120
dt=0.1
nt=0;
ntplot=1;

x=[1:nx]*dx;

% 2 phase slips creating wavelengts
eps=4
q=.8
% no phase slips
eps=4
q=1
% 2 phase slips deleting wavelengths
eps=4
q=1.23

y=.4*sin(1.23*x+eps*sin(x/(2*nwave)));

% random initial conditions
amp=.001
y=amp*randn(1,nx);
 
Y=y';
t=0

e=ones(nx,1);

D0=spdiags([e],0,nx,nx);
D2=1/(dx*dx)*(spdiags([e -2*e e],-1:1,nx,nx)+spdiags([e],nx-1,nx,nx)+spdiags([e],-nx+1,nx,nx));
D4=1/(dx*dx*dx*dx)*(spdiags([e -4*e 6*e -4*e e],-2:2,nx,nx)+spdiags([e -4*e],nx-2:nx-1,nx,nx)+spdiags([-4*e e],-nx+1:-nx+2,nx,nx));
Mt=1/dt*spdiags([e],0,nx,nx);

Dl=Mt+alpha*D4+2*alpha*D2-alpha*(R-1)*D0;
Dr=Mt-(1-alpha)*D4-2*(1-alpha)*D2+(1-alpha)*(R-1)*D0;

% transponiere mit ': yt=y'
delete(fig1)
fig1=figure(1);
lin1=line(x,y,'EraseMode','background');
axis([0 xmax -1 1]);
drawnow;

while t < tmax
Yq=(Y.*Y);
Yc=(Y.*Y.*Y);
rhs=Dr*Y-Yc;
Dlc=3*alpha*spdiags([Yq],0,nx,nx);
Dls=Dl+Dlc;
Y=Dls\rhs;
y=Y';

if rem(nt,ntplot) == 0
delete(lin1);
lin1=line(x,y,'EraseMode','background');
drawnow;
end

nt=nt+1;
t=t+dt;
end

