A fractional slope on the log log Bode plot has been observed in characterizing a certain type of physical phenomenon and has been called the fractal system or the fractional power pole. In order to represent and study its dynamical behavior, a singularity function method is presented which consists of cascaded branches of a number of pole-zero (negative real) pairs or simple RC section. The distribution spectrum of the system can also be easily calculated, and its accuracy depends on a prescribed error specified in the beginning. The method is then extended to a multiple-fractal system which consists of a number of fractional power poles. The result can be simulated by a combination of singularity functions, each representing a single-fractal system.< >