La ecuación de la difusión neutrónica describe la población de neutrones dentro de un reactor nuclear. Este trabajo está enfocado a reactores nucleares con geometría hexagonal. En primer lugar se estudia la ecuación de la difusión neutrónica estacionaria. Este es un problema diferencial de valores propios, llamado problema de los modos Lambda. Para resolver el problema de los modos Lambda se han comparado diferentes métodos en geometrías unidimensionales, resultando como el mejor, el método de elementos espectrales. Usando este método discretizamos los operadores en geometrías bidimensiones y tridimensionales, resolviendo el problema algebraico de valores propios resultante con el método de Arnoldi. La distribución de neutrones en estado estacionario se utiliza como condición inicial para la integración de la ecuación de la difusión neutrónica dependiente del tiempo. Se utiliza un método de Euler implícito para integrar en el tiempo. Cuando una barra de control está parcialmente insertada en un nodo aparece un comportamiento no físico de la solución, el efecto ``rod cusping'', que se corrige mediante la ponderación de las secciones eficaces con el flujo del paso de tiempo anterior. Para la resolución de los sistemas algebraicos que surgen en el método implícito, se utiliza un método de Krylov, y se evalúan diferentes estrategias de precondicionamiento. La primera consiste en el uso de la estructura a bloques definida por los grupos de energía para resolver el sistema. Además se proponen diferentes técnicas de aceleración para el esquema iterativo de bloques y un precondicionador utilizando esta estructura. Además se estudia un preacondicionador espectral, que hace uso de la información del subespacio de Krylov obtenida cuando se resuelve un sistema para precondicionar el siguiente sistema. También se proponen métodos exponenciales de segundo y cuarto orden para integrar la ecuación de difusión neutrónica dependiente del tiempo, donde la exponencial de la matriz del sistema tiene que ser multiplicada por un vector. Estos esquemas nos permiten trabajar sin construir explícitamente la matriz del sistema, y se comparan diferentes métodos para aproximar el producto de la exponencial matricial por un vector, tal como un método de Krylov, un método de Chebyshev y un método basado en los puntos de Leja. Surgen algunas situaciones en las que un conjunto de modos tiene que ser actualizado, como en los cálculos perturbados o en el uso de métodos modales. La actualización de estos de modos puede ser muy costosa cuando se utiliza el método de Arnoldi y, por esta razón, se proponen varios métodos basados en una iteracion de Newton, tales como el método de Newton a bloques modificado, y dos alternativas basadas en iterar con uno o dos subespacios. Como estrategia alternativa se propone un modelo de orden reducido para actualizar un conjunto de modos basados en la ``Proper Generalized Decomposition''. Este método obtiene una solución aproximada como una suma de funciones separables sobre todo el dominio, reduciendo el problema multidimensional a un conjunto de problemas unidimensionales.