Multicanonical ensemble sampling simulations have been performed to calculate the phase diagram of a Lennard-Jones fluid embedded in a fractal random matrix generated through diffusion- limited cluster aggregation. The study of the system at increasing size and constant porosity shows that the results are independent of the matrix realization but not of the size effects. A gas-liquid transition shifted with respect to bulk is found. On growing the size of the system on the high-density side of the gas-liquid coexistence curve there appears a second coexistence region between two liquid phases. These two phases are characterized by a different behaviour of the local density inside the interconnected porous structure at the same temperature and chemical potential.