For the mode matching calculation I was using the ETMY focal length that I found on Kiwamu's plot on the wiki page. 

Taking into account also the substrate, the focal length turns out to be

fl = ((n-1)*(1/R1 - 1/R2 + (n-1)d/(nR1R2)))^(-1) = -125.81 m

with n = 1.46071 (refraction index of fused silica at 532nm)

R1 = 5625 m (radius of curvature of the first surface)

R2 = 57.37 m (radius of curvature of the second surface)

d = 25mm (thickness)

The value of the focal length is sligthly different from the one I was using before in the calculation, but maybe it is enough to change the coupling.

The mode matching solution I found is very sensitive to the lenses position. 

The beam waist position can vary up to 20m varying by 1cm the first lens position, while it is slightly less sensitive to the second lens displacement.

As shown in the picture, along the green beam path there is also a 1m focal length lens. It's position is fixed, because it is along the IR transmetted beam path also. I tried to get a better solution without it, but I found that the waist position was still strongly dependent on one of the two lenses position, so it would not solve the problem to remove this lens.

I think that the main issue of this mode matching is related to the "space contraints", because the two lenses' positions can vary in a very small space, even though the green beam path on the table is quite long.

Eventually, I put the MM lenses found from this last simulation on the table, and it seems to work, since I've seen very strong 00 flashes. Unfortunately, while trying to maximize the alignment I broke it  and I have to do it again, but I feel confident!