A method is proposed to describe Fermi or Bose systems coupled to one or several heat baths composed of fermions and/or bosons. The method, called Coupled Equations of Motion method, properly includes non-Markovian effects. The approach is exact in the Full-Coupling approximation when only bosonic particles are present in the system and baths. The approach provides an approximate treatment when fermions are present either in the system and/or in one or several environments. The new approach has the advantage to properly respect the Pauli exclusion principle for fermions during the evolution. We illustrate the approach for the single Fermi or Bose two-level system coupled to one or two heat-baths assuming different types of quantum statistics (Fermion or Bosons) for them. The cases of Fermi system coupled to fermion or boson heat baths or a mixture of both are analyzed in details. With the future goal to treat Fermi systems formed of increasing number of two-level systems (Qubits), we discuss possible simplifications that could be made in the equations of motion and their limits of validity in terms of the system-baths coupling or of the initial heat baths temperatures.