Thread Subject: The unstableness of my total variation denoising

I am using the gradient descent method for total variation regularized reconstruction. Specifically, the following code is used for the time evolution:

     xi = xi + tau*( Num./ Den*lambda + fidelity);

where the Num./Den is the PDE flow computed using the upwind finite difference scheme and the lambda is the regularization parameter. Besides, the smoothing parameter for TV norm is set as: eps=1 and the tau=0.2;

I found that the solution becomes highly unstable, if I increased the lambda for strong regualrization. Personally, I think this is because the CFL condition is broken. However, I didn't see any literature mentioned this problem. Is there a mistake in the parameter selection such as tau and eps or is this a common phenomenon?

Another question is how to choose the eps? In my case, it seems that eps below one may make the Den too small and result in a very large flow at some pixels. And the whole image was ruinted. However, I remembered that a small eps will have better denoising result. So I don't know how to handle these two things simultaneously.
And I sincerely hope you guys can give me some suggestions.
Thank you so much in advance!

################################################
P.S.
Bellow is a full version of my code for the total variation reconstruction:

for kk = 1:n_iter
    tic

    % compute the fidelity term A'b - A'Ax
    q = compAx(xi,s,options);
    fidelity = b - q; % fidelity = A'b - A'Ax

        I_x = (xi(:,[2:nx nx])-xi(:,[1 1:nx-1]))/2;
        I_y = (xi([2:ny ny],:)-xi([1 1:ny-1],:))/2;
        I_xx = xi(:,[2:nx nx])+xi(:,[1 1:nx-1])-2*xi;
        I_yy = xi([2:ny ny],:)+xi([1 1:ny-1],:)-2*xi;
        Dp = xi([2:ny ny],[2:nx nx])+xi([1 1:ny-1],[1 1:nx-1]);
        Dm = xi([1 1:ny-1],[2:nx nx])+xi([2:ny ny],[1 1:nx-1]);
        I_xy = (Dp-Dm)/4;

       % compute flow
        Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
        Den = (ep2+I_x.^2+I_y.^2).^(3/2);

       %% evolve image by options.tau
        xi=xi+options.tau*(Num./Den*lambda + fidelity);

   toc
end