If there are only two players, you could use integer linear programming.
I'll introduce some zero-or-one integer variables to encode the (unknown) bidding functions (via a one-hot encoding). Let $x_{v,b}=1$ if player $0$'s bid is $b$ when their valuation is $v$, and $0$ otherwise; and $y_{w,c}=1$ if player $1$'s bid is $c$ when their valuation is $w$, and $0$ otherwise. In this way, if there are $n$ possible valuations, we obtain $2n^2$ zero-or-one integer variables. We will write down linear inequalities on these variables that characterize the optimality condition, then use an integer programming solver to find a solution that satisfies all of these inequalities.
Note that if player $0$'s valuation is $v$ and player $1$'s valuation is $w$, then player $0$'s payoff (in a first-price auction) is
$$\sum_{b>c} x_{v,b} y_{w,c} (v-b).$$
Therefore, if player $0$'s valuation is $v$, player $0$'s expected payoff is
$$\sum_w p_1(w) \sum_{b>c} x_{v,b} y_{w,c} (v-b),$$
where $p_1(w)$ represents the probability that player $1$ gets valuation $w$. Now, we need this to be at least what it would be if player $0$ used any other bid for valuation $v$, i.e., any other choice of $x_{v,\cdot}$'s. If there are $n$ possible valuations, there are only $n$ possible bids, so we obtain $n$ inequalities
$$\sum_w p_1(w) \sum_{b>c} x_{v,b} y_{w,c} (v-b) \ge \sum_w p_1(w) \sum_{b'>c} y_{w,c} (v-b'),$$
where we obtain one such inequality for each possible value of $b'$. (The sum on the left-hand-side is over $w,b,c$ that satisfy the condition $b>c$; the sum on the right-hand-side is over $w,c$ that satisfy the condition $c<b'$.)
This is a quadratic inequality. We'll turn it into a linear inequality using the techniques of Express boolean logic operations in zero-one integer linear programming (ILP). In particular, introduce new zero-or-one variables $t_{v,w,b,c}$, with the intention that $t_{v,w,b,c}=x_{v,b} y_{w,c}$. This can be enforced using the linear inequalities $t_{v,w,b,c} \ge x_{v,b} + y_{w,c} - 1$, $t_{v,w,b,c} \le x_{v,b}$, $t_{v,w,b,c} \le y_{w,c}$, $0 \le t_{v,w,b,c} \le 1$. Now the optimality condition becomes
$$\sum_w p_1(w) \sum_{b>c} t_{v,w,b,c} (v-b) \ge \sum_w p_1(w) \sum{b'>c} y_{w,c} (v-b'),$$
for each $v,b'$. We obtain one such inequality for each $v$ and each $b'$. These $n^2$ inequalities ensure that player $0$'s strategy will be optimal (given player $1$'s strategy).
Then, add similar inequalities to require player $1$'s strategy to be optimal (given player $0$'s strategy).
Also add the $n^2+n$ constraints that $0 \le x_{v,b} \le 1$ and $\sum_b x_{v,b}=1$ to require the $x$'s to correspond to a bidding function, and similarly for the $y$'s.
Finally, take all of these inequalities and feed them to an off-the-shelf integer linear programming solver. If it finds a solution, then you have found a Nash equilibrium for your game. ILP is NP-hard, so there is no guarantee that it will terminate in a reasonable amount of time, but in my experience often the solvers can handle surprisingly large systems of inequalities.
The same approach can handle all-pay auctions as well as first-price auctions.
In practice, you might be able to speed up the solver by adding some extra constraints to help it identify a solution faster. In particular, I have a hunch that without loss of generality the optimal bidding function can be assumed to be monotonic: if we increase my valuation, my bid should never decrease. We can add extra inequalities to the system to characterize this additional assumption, and this might narrow the search space and thus help the ILP solver find a solution more rapidly.
One big limitation is that this seems limited to two players (or a very small number of players). If there are multiple players, the number of variables and inequalities will explode, and this might not be an effective approach.
