I have found a much more efficient way of solving this problem on computer. Having asked the question I guess I should provide a brief account. However I feel like the algorithm is technical and not very enlightening, and so I'll be brief. Please leave comments if you would like more detail.

First we consider an important special case: Is there an efficient algorithm to find all subexpressions of a fixed (reduced) expression $\underline{w} = [s_1, \dots, s_m]$ whose product is the identity? In the language of the question, how does one calculate $\Pi(\underline{w}, id)$?

Of course, there is one canonical such subexpression: $(id,id,…,id)$. Moreover, given any subexpression for the identity we can apply "cancellation moves" to get to this expression.

By a "cancellation move" I mean the following: write $\pi[i,j]$ for the product $t_it_{i+1}…t_j$. Suppose that

$s_i \pi[i+1,j-1] = \pi[i+1,j-1] s_j$

then a "cancellation move" is one of the following:

$\pi = [\dots ,s_i, \dots, id, \dots] \mapsto \pi ' = [\dots, id,…,s_j, \dots]$

$\pi = [\dots ,s_i, \dots, s_j, \dots ]\mapsto \pi' = [\dots ,id, \dots, id, \dots]$

where $\pi$ remains unchanged except at the $i^{th}$ and $j^{th}$ place. Obviously, applying a cancellation move to a subexpression does not change the product. Also, cancellation moves either increase the number of $id$'s or increase the number of $id$'s to the left. It follows that by repeatedly applying cancellation moves to any subexpression for the identity we end up at the canonical subexpression.

For example, consider the Coxeter group with simple reflections $s$, $t$ and $u$ and only braid relation $(st)^3 = id$ (so that $su$ and $tu$ have infinite order). Let $\underline{w} = [s,t,s,u,t,s,t]$. Then the following is a sequence of cancellation moves:

$[s,t,s,id,t,s,t] \mapsto [id,t,s,id,id,s,t] \mapsto [id,t,id,id,id,id,t] \mapsto [id, id, id, id, id, id, id]$.

Now one can reverse this, and define "reverse cancellation moves". Applying all possible reverse cancellation moves to the canonical subexpression yields all subexpressions for the identity. (This is easily programmed.)

We now turn to the general case. The above "cancellation moves" make sense for any subexpression and one can show that again there is a canonical subexpression. (It is characterised by the fact that the subexpression is reduced, and the $t_i \ne id$ occur as far to the right as possible.) Again, once one has this canonical subexpression then it is a simple matter to reconstruct all subexpressions by reverse cancellation moves.

Of course this begs the question: can one find this canonical subexpression efficiently? Yes. Firstly, consider $R(x)$ the right descent set of $x$, and choose $s \in R(x)$ such that s occurs as far as possible to the right in $\underline{w}$. Now, delete everything to the right and including the right-most occurrence of $s$ and repeat with $x$ replaced by $xs$. (This also gives a reasonably efficient algorithm to decide whether $x \le y$ in the Bruhat order.)

(One can ask what this canonical expression "means". Here is one possible explanation: because Schubert varieties are normal the fibre of the Bott-Samelson resolution over point corresponding to $x \le w$ is connected and hence its homology is one-dimensional in degree zero. Now subexpressions also index BB-cells, and hence a basis for the cohomology of the fibre. This "canonical subexpression" corresponds to a generator of the one-dimensional degree zero part.)

For example, take $W = S_7$ (with simple reflections = simple transpositions) and take the following reduced expression for the longest element:

$w_0 = 121321432154321654321$

then there are 6408 subexpressions with product the identity. On my laptop the above algorithm takes 3.33 seconds to find them. A brute force attack takes 93 seconds.

Here is another example which I actually cared about in my calculations. Take the reduced expression

$\underline{w} = 13572613574352461357$

in $S_8$. There are 80 subexpressions with product the longest element in the standard parabolic subgroup generated by all simple transpositions except 4. The above algorithm takes 0.03 seconds to find all of them, whereas the naive approach takes roughly 15 seconds. (So here the above algorithm is 500 times as fast.)

I imagine that the differences become (even) more pronounced with longer words.