#
# Numerische Lsg. der KdV Gleichung
# 
# u_t + ep u u_x + mu u_xxx = 0
#
# Alex Weisse 9/6/08
#

# Parameter der KdV Gl.
ep = -6.0;
mu = 1.0;

dx = 0.2;
dt = 0.001;

# Gittergroesse
x = -10:dx:10;
t = 0:dt:4;

nx = length(x);
nt = length(t);

# Anfangsbedingung
u = zeros(nt, nx);
u(1,:) = -6*sech(x-5).^2;

# erster Zeitschritt
u(2,:) = u(1,:)\
   -ep/6*dt/dx*(shift(u(1,:),1)+u(1,:)+shift(u(1,:),-1))\
        .*(shift(u(1,:),1)-shift(u(1,:),-1))\
   -mu/2*dt/(dx^3)*(shift(u(1,:),2)-2*shift(u(1,:),1)+2*shift(u(1,:),-1)
        -shift(u(1,:),-2));

# weitere Zeitschritte
for j = 2:nt-1
  u(j+1,:) = u(j-1,:)\
     -ep/3*dt/dx*(shift(u(j,:),1)+u(j,:)+shift(u(j,:),-1))\
         .*(shift(u(j,:),1)-shift(u(j,:),-1))\
     -mu*dt/(dx^3)*(shift(u(j,:),2)-2*shift(u(j,:),1)+2*shift(u(j,:),-1)\
         -shift(u(j,:),-2));
endfor

# graph. Darstellung
vtk_clear;
[cx,ct] = meshgrid(x,t);
vtk_ylabel("10 t");
vtk_zlabel("u(x,t)");
vtk_surf(cx,10*ct,u);
