As noted in the other answer, your problem is not using both outputs from gradient. However, here I want to note important issues with the formulation.
First, your discrete problem can be written as a standard linear least squares problem. That is, your objective function is (half of)
$$E[u] = \|u-f\|^2 + \lambda(\|\Delta_xu\|^2+\|\Delta_yu\|^2)$$
where now $u$ and $f$ are vectors created by flattening the corresponding image arrays (i.e. u=U(:) and f=F(:)).
This is a sum of squares, which is equivalent to the least-squares residual
$$
\begin{bmatrix} I \\ \sqrt{\lambda}D \end{bmatrix}
u \approx \begin{bmatrix} f \\ 0 \end{bmatrix}
\,,\,
D = \begin{bmatrix} D_x \\ D_y \end{bmatrix}
$$
where $I$ is the identity matrix, $D$ is a discrete gradient operator.
The least squares solution of the system $Au=b$ can then be computed as u=A\b, using Matlab's backslash operator.
Second, for this problem the discrete gradient would typically be computed using forward differences (rather than the centered differences used by gradient). This is because the least squares solution of the above problem is
$$(I+\lambda{D^T}D)u=f$$
With 1st order differences, the regularization matrix $D^TD$ is then a 2nd order discretization of the Laplacian $\nabla^2u$. The equation above is then a discrete Poisson equation, and moreover is consistent with the continuous (variational) form of the optimization problem.
If all of the above is done with sparse matrices, then the solution is very efficient, and with no iteration required. As a side note, this can all be done in Matlab using gridfit (with the "springs" option).
For completeness, below is an example Matlab implementation.
[ny,nx]=size(F);
Dfnc=@(sz)spdiags(ones(sz(1),1)*[-1,1],0:1,sz(1)-1,sz(1));
Dy=kron(eye(nx),Dfnc(ny,nx)); Dx=kron(Dfnc(nx,ny),eye(ny)); D=[Dx;Dy];
% Fx=diff(F,[],2); Fy=diff(F,[],1); Derr=norm(D*F(:)-[Fx(:);Fy(:)]) % Derr=0
I=speye(ny*nx); A=[I;lam*D]; b=[F(:);zeros(size(D,1),1)];
U=zeros(ny,nx); U(:)=A\b;
Best Answer
I don't know algorithms off the top of my head, but I would start by having a look at Survey of image registration techniques