A BFS starting from the solved cube is the way to go. However, you don't use the complete cube but instead ignore all the aspects that are not affected by the current stage of the algorithm.
In the first stage, where the edge orientations get solved, the only aspect of the cube that you look at are those edge orientations. Ignore the corners, and think of all the edges as being identical pieces so that their permutation is irrelevant. Of course the edges still move around when a move is done, but you cannot distinguish one edge piece from another. Only their orientation matters, so you can still see whether an edge is flipped or not.
Do a BFS with such a simplified cube to build your pruning table. The state of the cube can be uniquely encoded into an 11-bit integer, where each bit is the orientation of an edge at some location on the cube. That integer is the index into your pruning table. The table will store the distance from solved for each state of this simplified cube. Each time the BFS encounters a cube state it has not seen before (the associated table entry has not yet been filled), it then stores the current distance from start into the table. It is basically generating a God's Algorithm table for this simplified cube.
You will need to define what the edge orientation is on your cube. It depends on which pair of opposite faces gets restricted in $G_1$. Let's say it is F and B. On the simplified edge-orientation-only cube any moves in G1=<U,D,R,L,F2,B2> do not flip any pieces (though they are still moved around). The moves F, B, F', and B' however will each not only move four edges around, they also flip those four moving edges. With this definition you can build your pruning table through a BFS search.
When you want to solve a normal cube, and use this pruning table to solve the first stage in the algorithm, you will need to be able to look up the normal cube in your pruning table. So you need to extract the edge-orientation-index from it. For any edge, imagine bringing it to its home location using only moves from $G_1$. Once it is in its home location you can easily tell if it is flipped or not by whether its colours match the face centres. Since the $G_1$ moves did not affect its orientation, you now know whether the edge is flipped even if it was not at its home location.
When you have determined the orientation of all the edges of your cube, you can set the simplified edge-orientation cube to that state, use the pruning table to find a solution for this stage, and then apply this solution to the original cube.
The other stages of the algorithm can be done in essentially the same way, generating one or more pruning tables using a simplified version of the cube that only encodes the aspects being solved in that stage. The most difficult aspect to encode into a pruning table index will be what happens to the corner permutation in the third stage.
For each stage you can either generate a single pruning table, so that solving that stage for a real cube is essentially a God's Algorithm set of table lookups, or you can use several pruning tables, in which case solving that stage will involve a search using those tables (i.e. IDA*, iterative deepening A* search).
I presume you have found the page on my website about the Thistlethwaite algorithm. Another page that might be helpful is the page about computer cubing which has brief descriptions of the techniques involved.
Best Answer
not sure who told you that but according to https://cube20.org/ there are 531,653,418,284,628 positions solvable in exactly 13 moves (and no less).