To evaluate the availability and performance of pencil beam algorithm for fast neutrons radiotherapy dose calculations and to achieve improvements in universality, speed and inhomogeneities corrections. Pencil beam (PB) method uses integration of dose kernel cross field size to estimate dose distribution. In order to calculate dose kernels we used general purpose Monte-Carlo (MC) code MCNP. The same program was used also for PB benchmarking. Multi-group approach used for kernels calculation virtually allows calculation for arbitrary spectrum. Nevertheless we took four common spectra to investigate accuracy and to establish according corrections. Several approaches for inhomogeneities corrections have been investigated. The output of MC calculations is discrete kernel data, which does not allow fast calculations. To overcome this limitation we investigated different dose kernels approximations and integrations. PB dose kernel library was established. Comparison of PB calculation with reference algorithm shows adequate precision for arbitrary neutron spectra with inhomogeneities correction, while calculations are fast enough for optimization tasks. Secondary photon dose distributions comparison give worse still promising results. Researches of this paper show benefits and limitations of PB method for neutron teletherapy treatment planning. Calculations made during research allow saying that results of PB method are generally acceptable. Still additional data are needed to evaluate our results, probably voxel phantom MC calculations and/or beam data measurements.