Hello Annisa,
the function currently only implements the overlap of two species.
If you need more maybe you can run activityDensity on all relevant species, save the output (which it time in radians, between 0 and 2pi). Then use that as input to overlap::densityPlot(), which returns a data frame with x and y coordinates for further plotting.
You could then combine these data frames for multiple species and plot with ggplot2, for example.
This is just off the top of my head and untested. It also doesn't give overlap coefficients. If you need those you could calculate them with overlap::overlapEst(), but that would still be pairwise only, not for all species at once
Feel free to let us know how it's going.
Best regards,
Jürgen