Since they consider $A'$ and $\lambda'$ to be local functions, they can only be those constructed form $A,\partial A, \lambda, \partial \lambda$ and possibly higher order derivatives. Since there's only the constraint on the correspondences on gauge transformation, they probably guessed some forms of combinations of those variables, and find one possible solutions. They did not claim the map is unique, so you are probably right in that the single equation system is under-constrained to completely determine the map.
By treating $A$ and $\partial A$ as independent variables, and guessing $A'_i(A)=f_{ij}(A) A_j$ and $\lambda'(A,\lambda)=g_i(A,\lambda)A_i$. Setting $\lambda=0$ and $A=0$, the constraint becomes $f_{ij}(A)\partial_j \lambda-\partial_i (g_j(A,\lambda) A_j)=-\frac12\theta^{kl}\{\partial_k\lambda,\partial_l A_i\}.$ If we further guess $g_j(A,\lambda)A_j=g_{ij}\partial_i\lambda A_j$, replacing $\partial_j\lambda$ by $A_j$ leads to something mixture of $\theta^{kl}$ and $g_{kl}$ in the expression of $A'$ which now depends on $A$ and $\partial A$.
Some other choices on setting $\lambda$, $A$ and $\partial \lambda$, $\partial A$ can then be used to simultaneous determine $g_{kl}$ and $f_{ij}(A)$. So my guess is that they did some trial-and-error and figured out one (probably non-unique) map.