function KriegD

    globals;

    ux = (u(2:N+1)-u(1:N))/h;
    Di = 1./sqrt(ux.^2 + ep^2);
    D  = zeros(N+1);

    for i=1:N
        D(i,i+1) = -Di(i);
        D(i+1,i) = -Di(i);
    end

    for j=1:N-1
        D(j+1,j+1) = Di(j) + Di(j+1);
    end

    D(1,1)      = Di(1);
    D(N+1,N+1)  = Di(N);

% ohne Schleifen:
%   D = diag([Di;0]) ...
%     + diag([0;Di]) ...
%     - diag(Di,-1)  ...
%     - diag(Di,+1);

end