function u_kp1 = debruitage(b, u_k, lambda, Dx, Dy, epsilon)

    N = length(b);
    u_kp1 = zeros(size(u_k));

    for i=1:3

        coeffs = 1 ./ sqrt( (Dx * u_k(:,:,i)).^2 + (Dy * u_k(:,:,i)).^2 + epsilon );
    
        W = spdiags(coeffs, 0, N, N);
    
        A = speye(N) + lambda * (Dx' * W * Dx + Dy' * W * Dy);
    
        u_kp1(:,:,i) = A \ b(:,:,i);

    end

end