The long-term evolution of the internal, strong magnetic fields of neutron stars needs a specific numerical modelling. The diversity of the observed phenomenology of neutron stars indicates that their magnetic topology is rather complex and three-dimensional simulations are required, for example, to explain the observed bursting mechanisms and the creation of surface hotspots. We present MATINS, a new three dimensions numerical code for magneto-thermal evolution in neutron stars, based on a finite-volume scheme that employs the cubed-sphere system of coordinates. In this first work, we focus on the crustal magnetic evolution, with the inclusion of realistic calculations for the neutron star structure, composition and electrical conductivity assuming a simple temperature evolution profile. MATINS follows the evolution of strong fields (1e14-1e15 Gauss) with complex non-axisymmetric topologies and dominant Hall-drift terms, and it is suitable for handling sharp current sheets. After introducing the technical description of our approach and some tests, we present long-term simulations of the non-linear field evolution in realistic neutron star crusts. The results show how the non-axisymmetric Hall cascade redistributes the energy over different spatial scales. Following the exploration of different initial topologies, we conclude that during a few tens of kyr, an equipartition of energy between the poloidal and toroidal components happens at small-scales. However, the magnetic field keeps a strong memory of the initial large-scales, which are much harder to be restructured or created. This indicates that large-scale configuration attained during the neutron star formation is crucial to determine the field topology at any evolution stage.