Computing pure strategy Nash equilibria in finite games I am trying to compute the (pure strategy) Nash equilibria of some discrete auctions.

More precisely, let us define the strategy of each player as a function mapping from every valuation that they might have to their bid (i.e. the 'bidding function'). Let us suppose that each player's valuation is drawn from some finite set and that their bid must belong to this same (finite) set. I am interested in finding the set of bidding functions, one for each player, such that each player's bidding function is optimal given the bidding functions of all the other players.

If it makes things easier, we can assume that valuations are symmetric (i.e. each player's valuation is generated by the same probability distribution) and that there are only two players. Ideally, however, we would proceed without these simplifications. I am interested in computing the equilibrium for the first price sealed bid and all pay auctions (in the latter, you pay your bid even if you lose; in the former, you don't.)

I have considered writing down the discrete auction game in normal form and finding the equilibria using software like Gambit. However, this would seem tricky since the strategy space is so large.  For example, if a player chooses bids from $ \{1,...,10\} $ and draws values from $ \{1,...,10\} $, then already they have $10^{10}$ pure strategies.

Does anyone have any ideas about how to proceed here?
If $10^{10}$ pure strategies is too large for Gambit, it'll likely be too large for any other software as well. In the comments, I suggested that you could first compute the expected payoffs for each player in a separate program such as Matlab/R/Excel, and then export the resulting matrix to Gambit where the BNEs can be calculated. This way, you may be able to save some resources for Gambit on converting the game into its Bayesian normal form. 

If you're okay with significantly reducing the dimension of your problem, then you can rely solely on Gambit (or, Game Theory Explorer (GTE), which is now a part of the Gambit project). In the following example done in GTE, I assumed that bidder $i$'s valuation and strategy spaces are binary: $V_i=S_i=\{0,1\}$ for $i=1,2$. The values are assumed to be uniformly distributed. 

First, you enter the game in its extensive form. In the figure below, the four branches at the top correspond to Nature's choices that determine the valuations. From left to right, the branches represent the valuation profiles $(v_1,v_2)=(0,0)$, $(0,1)$, $(1,0)$, $(1,1)$ respectively, each occurring with probability $1/4$. 



After you finish entering the extensive form game, press the "Matrix Layout" tab to see the (automatically calculated) Bayesian normal form. There's a gear button at the top that will compute the BNE of the game.



The result looks like the following. In this case, the result was returned in seconds (in part thanks to the relatively small dimension of the problem). The unique BNE, in which each player bids zero regardless of their value, was found.  



You can of course try to increase the size of the value/strategy space. This method would also allow you to specify asymmetric distributions of bidder values. However, GTE only solves two-player games. But Gambit seems to be able to solve games with more than two players.