Raising the order of the multipole expansion is a feasible approach for improving the accuracy of the treecode algorithm. However, a uniform order for the expansion would result in the inefficiency of the implementation, especially when the kernel function is singular. In this paper, a $p$-adaptive treecode algorithm is designed to resolve the efficiency issue for problems defined on a general domain. Such a $p$-adaptive method is realized through i). conducting a systematical error analysis for the treecode algorithm, ii). designing a strategy for a non-uniform distribution of the order of multipole expansion towards a given error tolerance, and iii). employing a hierarchy geometry tree structure for coding the algorithm. The proposed $p$-adaptive treecode algorithm is validated by a number of numerical experiments, from which the desired performance is observed successfully, i.e., the computational complexity is reduced dramatically compared with the uniform order case, making our algorithm a competitive one for bottleneck problems such as the demagnetizing field calculation in computational micromagnetics.