%initial info

n = 100;
dx = 1/n;
dt = 0.00001;
T = 0.05;
x = dx:dx:1-dx;

v0 = 10^6*x'.^5.*(1-x)'.^20; % initial velocity

c = 10;
d = 1;

K = toeplitz(sparse([2 -1 zeros(1,n-3)]));
S = toeplitz(sparse([6 -4 1 zeros(1,n-4)]));

r = c*dt/dx;
R = d*dt/(dx^2);

% Leapfrog

ULF = zeros(n-1, round(T/dt));
ULF(:,2) = ULF(:,1) + dt*v0;
M = 2*speye(n-1) - r^2*K - R^2*S;

for m=3:round(T/dt)
	ULF(:,m) = M*ULF(:,m-1) - ULF(:,m-2);
end

% One can use explicit first-order methods like Euler or Runge-Kutta on 
% the pair (u,u'), but stability would require dt/(dx)^4 < 1.