How to Generate Münchhausen Numbers in High Radices?
Meet-in-the-middle

One approach is to use a meet-in-the-middle algorithm.  This will reduce the running time from  $2b^b$ to something in the ballpark of $\sqrt{2b^b}$.

In particular, let your candidate number be
$$A = a_b b^b + \dots + a_1 b + a_0.$$
We'll write this in the form $A=A^* + A^\dagger$ where
$$\begin{align*}
A^* &amp;= a_b b^b + \dots + a_{\lfloor b/2 \rfloor+1} b^{\lfloor b/2 \rfloor+1}\\
A^\dagger &amp;=a_{\lfloor b/2 \rfloor} b^{\lfloor b/2 \rfloor} + \dots + a_0.
\end{align*}$$
Finally, let $s(a_i b^i + \dots + a_0) = a_i^{a_i} + \dots + a_0^{a_0}$.  We'll look for all solutions to the equation
$$A^* + A^\dagger = s(A^*) + s(A^\dagger).$$
This is equivalent to
$$A^* - s(A^*) = s(A^\dagger) - A^\dagger.$$
We'll compute a table of all possible values for $s(A^\dagger)-A^\dagger$, where $A^\dagger$ varies over all possibilities of the form above.  Then, we'll compute $A^*-s(A^*)$, for each $A^*$ of the form above, and look it up in the table to check for a match.  Each match yields a solution to the equation above, and thus a valid Münchhausen number.  This can be used to enumerate all Münchhausens.

For $b=10$, the running time becomes something like $10^5$ steps of computation, which should be no problem.  However, this is entirely infeasible for $b=35$.  The algorithm is still exponential-time.

Linear programming

Another possible idea is to use integer linear programming.

Let $x_{i,d}$ be a 0-or-1 variable, with the intended meaning that $x_{i,d}=1$ means that $a_i = d$.  Then we can obtain a bunch of linear inequalities that characterize a correct solution.  In particular, we have
$$\sum_{i,d} x_{i,d} d b^i = \sum_{i,d} x_{i,d} d^d.$$
We also have $\sum_d x_{i,d} = 1$ for each $i$.

Next, given a lower and upper bound $L,U$, we could use integer linear programming to search for a valid solution in the range $[L,U]$, by adding the inequality
$L \le \sum_{i,d} x_{i,d} d b^i \le U$ and looking for a valid solution.  To enumerate all solutions, start with a broad range $[0,2b^b]$, then when you find a valid solution, split the current range in two and recurse.  Then, use an off-the-shelf ILP solver to find solutions.

I don't expect this to be very effective, but you could augment it with additional inequalities that might help the ILP solver run faster.  In particular, for each prime number $p$ less than $2b^b$, we have
$$\sum_{i,d} x_{i,d} d b^i \equiv \sum_{i,d} x_{i,d} d^d \pmod p$$
This can be encoded via the linear equality
$$\sum_{i,d} x_{i,d} d b^i = \sum_{i,d} x_{i,d} x_{i,d} d^d + y_p p$$
where $y_p$ is an integer variable.  Note that $d b^i$ and $d^d$ are integer constants, so you can reduce them modulo $p$. Thus we get the linear equality
$$\sum_{i,d} \alpha_{i,d,p} x_{i,d} = \sum_{i,d} \beta_{d,p} x_{i,d} + y_p p$$
where $\alpha_{i,d,p} = d b^i \bmod p$ and $\beta_{d,p} = d^d \bmod p$ are constants you can precompute.  You can do the same if $p$ is a power of a prime.

So, you might add linear inequalities of this form for several well-chosen prime powers (say, all of the prime powers below some bound; or all small prime numbers, plus the prime power factors of $b$, $b-1$, $b+1$, $b^2$, $b^2-1$, $b^2+1$).  One could hope that these help the linear solver find solutions faster.

I don't know whether this will be effective, but it is something you could try.  It would not be too hard to implement and see how well it works.