Surface acoustic waves (SAWs) enjoy profound importance across many disciplines, and one of their most prominent applications are the vast ranges of SAW devices made of piezoelectric multilayers. Thus a stable and efficient algorithm of analysing key SAW parameters in arbitrary layered media is of wide interest. This paper introduces such an algorithm based on the spectral collocation method (SCM). It firstly explains the fundamental governing equations and their numerical calculations via the SCM in a self-contained way, and then demonstrates the technique on the well-studied ZnO/diamond/Si material system, where the key factors of phase velocity, electromechanical coupling coefficient and effective permittivity are evaluated and discussed in detail. Compared to the widely employed root-finding approach, it is shown that the SCM is intuitive to formulate and code, and is highly accurate; it delivers all modes at once, and does not suffer from numerical instabilities. The establishment of the method on this simple example also implies potential applications to more general material and geometry types that could be difficult for the conventional approaches.