    N  = 20;
    om2 = 1;
    c  = 1;
    T  = 10;
    K  = 1000;
    ta = T/K;

    N1 = N+1;
    N2 = N+2;
    h  = 1/N1;

    D  = spdiags(kron([-1,1],ones(N1,1)),[-1,0],N1,N1);
    D  = D/h;

    g  = [zeros(N1,1),zeros(N1,1),-ones(N1,1)];

    s  = linspace(0,1,N2);
    u1 = s(2:end)';
    u2 = zeros(N1,1);
    u3 = zeros(N1,1);
    u  = [u1,u2,u3];

    ut = zeros(N1,3);
    U = [u;ut];

    for k=1:K
        u  = U(1:N1,:);
            u1 = [0;u(:,1)]; u2 = [0;u(:,2)]; u3 = [0;u(:,3)];
            plot3(u1,u2,u3); grid on; axis([-1 1 -1 1 -2 0]); drawnow;
        du = sqrt(sum((D*u).^2,2));
        a  = om2*(1-1./du);
        A  = D'*spdiags(a(:),0,N1,N1)*D;
        B  = [sparse(N1,N1),speye(N1);-A,-c*speye(N1)];
        U  = (speye(2*N1) - ta*B) \ (U + ta*[zeros(N1,3);g]);
    end
