Accurate modelling of wave-structure interaction is essential for a successful and reliable design of offshore structures and their ancillary systems. However, the fidelity of mathematical models is challenged when nonlinearities become significant, i.e. when the floater operates in severe sea states and/or it shows large dynamic responses compared to the incoming wave. This is normally the case for wave energy converters (WECs) because large motions are desirable for better power extraction, while conventional offshore structures are usually designed to prevent large responses and ensure stability. Fast computation is a further mandatory requirement for mathematical models that are used for design purposes. This is particularly true for wave energy applications, since extensive parametric studies are needed to minimize the cost of energy, i.e. increase power extraction capabilities while limiting capital and operational expenditures. Furthermore, model-based control strategies, essential for substantial increase of the WEC performance, require representative mathematical models. The requirements of accuracy and low computational time are usually contrasting, so an appropriate compromise should be defined. This paper proposes a nonlinear Froude-Krylov (NLFK) force model for axisymmetric floaters. The symmetry of the geometry (common for WECs) is exploited to achieve a low computational time (about real-time computation). Furthermore, since NLFK forces are the main nonlinearities for such devices, the obtained accuracy is higher than linear models. An open-source Matlab demonstration toolbox is introduced, called NLFK4ALL (doi: 10.5281/zenodo.3544848), which provides a ready-to-use implementation of the NLFK approach in six degrees of freedom. Accuracy and computational aspects of the method are here discussed, such as the complexity of the analytical description of the intersection between free surface elevation and the floater, and the impact of tolerances on the convergence of the numerical integration.