A new algorithm for direct numerical computation of aberration coefficients characterizing the imaging properties of charged particle optical devices is presented. The algorithm is based on numerical integraton of a system of differential equations derived from the equations of motion for a single particle. The proposed algorithm is free from a priori assumptions about electrostatic/magnetostatic field symmetries and is useful in the calculation of aberration coefficients around an arbitrary optic axis. The algorithm is illustrated by a calculation of aberration coefficients for a simple electrostatic lens.