from astropy.io import fits
from astropy.table import Table

my_table1 = Table(data=dict(a=[1, 2, 3]))
my_table2 = Table(data=dict(b=[4, 5, 6]))
hdu_list = fits.HDUList([
    fits.PrimaryHDU(), # need an empty primary HDU
    fits.table_to_hdu(my_table1),
    fits.table_to_hdu(my_table2),
])
hdu_list.writeto('tables.fits')