We study a family of $H^m$-conforming piecewise polynomials based on the artificial neural network, referred to as the finite neuron method (FNM), for numerical solution of $2m$-th-order partial differential equations in $\mathbb{R}^d$ for any $m,d≥1$ and then provide convergence analysis for this method. Given a general domain Ω$⊂\mathbb{R}^d$ and a partition $\mathcal{T}_h$ of Ω, it is still an open problem in general how to construct a conforming finite element subspace of $H^m$(Ω) that has adequate approximation properties. By using techniques from artificial neural networks, we construct a family of $H^m$-conforming functions consisting of piecewise polynomials of degree $k$ for any $k≥m$ and we further obtain the error estimate when they are applied to solve the elliptic boundary value problem of any order in any dimension. For example, the error estimates that $‖u−u_N‖_{H^m(\rm{Ω})}=\mathcal{O}(N^{−\frac{1}{2}−\frac{1}{d}})$ is obtained for the error between the exact solution $u$ and the finite neuron approximation $u_N$. A discussion is also provided on the difference and relationship between the finite neuron method and finite element methods (FEM). For example, for the finite neuron method, the underlying finite element grids are not given a priori and the discrete solution can be obtained by only solving a non-linear and non-convex optimization problem. Despite the many desirable theoretical properties of the finite neuron method analyzed in the paper, its practical value requires further investigation as the aforementioned underlying non-linear and non-convex optimization problem can be expensive and challenging to solve. For completeness and the convenience of the reader, some basic known results and their proofs are introduced.