I want to minimize the functional of teh Blind Deconvolution model as given in: Total Variation Blind Deconvolution by Chan and Wong.

Their model is given by:

$$ z = h \ast u + \eta $$

Where $ \ast $ is the convolution operator, $ h $ is the blurring kernel, $ u $ is the sharp noiseless image and $ \eta $ is additive white gaussian noise (AWGN).

The functional to be minimized is given by (Assuming the Blurring Kernel is known):

$$ \underset{u}{\min} f \left( u \right ) = \underset{u}{\min} \left \{ \frac{1}{2} {\left \| h \ast u - z \right \|}^{2}_{L_{2}} + \alpha \int \left | \nabla u \right | dx dy \right \} $$

Where $ \alpha $ is the smoothing term.

Yet, in Blind Deconvolution the Kernel isn't known and the minimization functional is given by:

$$ \underset{u}{\min} f \left( u \right ) = \underset{u}{\min} \left \{ \frac{1}{2} {\left \| h \ast u - z \right \|}^{2}_{L_{2}} + {\alpha}_{1} \int \left | \nabla u \right | dx dy + {\alpha}_{2} \int \left | \nabla h \right | dx dy \right \} $$

Now, The Euler Lagrange equations I calculated and given in the article are:

$$\begin{align} \frac{\delta L}{\delta h} & = \left ( u \ast h - z \right ) \ast u \left( -x, -y \right ) - {\alpha}_{2} \nabla \cdot \left( \frac{\nabla h}{\left | \nabla h \right |} \right) \\ \frac{\delta L}{\delta u} & = \left ( u \ast h - z \right ) \ast h \left( -x, -y \right ) - {\alpha}_{1} \nabla \cdot \left( \frac{\nabla u}{\left | \nabla u \right |} \right) \end{align}$$

What I'm not sure about is how can I use it to solve the problem.
At the article they suggest the alternating method, namely once solve for $ h $ and then for $ u $.

Yet I don't see how to write in in MATLAB code (Or any other pseudo code).
It should be some kind of a Gradient Descent step.