Instead of using the gradients/derivative, the update equation can be derived using the KL divengence.
Given the equation (21) in the question, the ELBO can be written as a function of $q_j$ without expanding the integration.
$ ELBO(q_j) = E_j[E_{-j}[\text{log}p(z_j,z_{-j},x)]] - E_j[\text{log}q_j(z_j)] + \text{constant} $
Where the first term is derived by writting the first term in (21) using the iterated expectation. The second term is derived by decomposing the second term in (21) as $q_j$ are independent with each other.
Actually, it's the same as the equation (22) in the question, just in another form.
Then, it can be found that $E_j[E_{-j}[\text{log}p(z_j,z_{-j},x)]] - E_j[\text{log}q_j(z_j)]$ is the negative KL divengence between $q_j$ and $\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}$.
So, in order to maximize the ELBO, we have to minimize the KL divengence between $q_j$ and $\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}$. The minimal happens when $q_j$ is equal to the $\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}$.
So, we can update the $q_j$ with:
$q^*_j\propto{\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}}$
The full distribution of $q^*_j$ can be get by normalizing according to the proportional relationship.
$q^*_j = \frac{ \text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\} }{ \int_{z_j}\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\} }$
But in reality, in order to update the variational parameter, the full distribution is usually not necessary to calculate if the conjugate exponential family distribution is used.
The core idea is explained above. For more detail and extension about how the exponential family is used in VI, refer the following link: