The recent discovery of neutron stars (NSs) with current gravitational-wave experiments gives us an unprecedented opportunity to probe the dynamical behaviour of compact stars. The viscous (and thermal) dissipation plays a fundamental role in the dynamical equations of motion of NSs. It has a damping influence on the oscillation modes and directly determines the minimum period of pulsars that is expected to be observed. We have studied radial pulsations of plausible NS models with various families of EoS in the presence of viscosity (and thermal conductivity) in the neutron matter. We have shown that the stellar pulsation equations can be cast in a nearly Sturm–Liouville form (just like the non-dissipative systems) and converted to a system of finite difference equations for numerical evaluation. We chose a second-order accurate differencing scheme so the resulting system of equations emerges as a tridiagonal matrix eigenvalue problem. Since radial oscillations do not couple to gravitational radiation, it is relatively easy to numerically solve the eigenvalue problem that leads to a discrete set of oscillation frequencies. In the absence of any dissipative process, the oscillation spectrum of a stable stellar model forms a complete set. It is therefore possible to describe any arbitrary periodic radial motion of a NS as a superposition of its various eigenmodes.