A semi-implicit multi-layer spherical spectral method for simulating stellar core convection is described. The fully compressible three-dimensional hydrodynamic equations with rotation and energy generation are solved. Prognostic variables are expressed as finite sums of spherical harmonics in the horizontal directions and handled by the finite difference method in the radial direction. The stratified approximation is used to simplify the nonlinearity to quadratic. A multi-layer scheme is employed to overcome the time step problem arising from shrinking grid sizes in the physical space near the center of the star. Despite of the different spectral truncations in different layers, round-off conservation of the total mass and total angular momentum of the whole domain can be maintained, and were confirmed numerically. The code is parallelized; with 12 processors the speedup factor is about 9. The solutions of model core convection with and without rotation are discussed.