The neutron diffusion equation describes the neutron population in a nuclear reactor core. This work deals with this model for nuclear reactors with hexagonal geometries. First, the stationary neutron diffusion equation is studied. This is a differential eigenvalue problem, called Lambda modes problem. To solve the Lambda modes problem, different methods have been compared in one-dimensional geometries, resulting as the best one the Spectral Element Method. The operators are then discretized using this scheme in two- and three-dimensional geometries, and the resulting algebraic eigenvalue problem is solved with the implicit restarted Arnoldi method.

Once the solution for the steady state neutron distribution is obtained, it is used as initial condition for the time integration neutron diffusion equation. Initially, a one step backard Euler method is used to discretize this equation in time. The transients to test the behaviour of the method are based on moving the control rods of the reactor, simulating an accident where a control rod is ejected and a scram is intialized to control the power evolution. An unphysical behaviour appears when a node is partially rodded, the rod cusping effect, which is corrected by weighting the cross sections with the flux of the previous time step. Good results are obtained for the power evolution and the rod cusping effect is corrected for different transients. To solve the algebraic systems arising in the backward method, a Krylov method is used, and different preconditioning strategies are tested. The first one consists of using the block structure obtained by the energy groups to solve the system by blocks, and different acceleration techniques for the block iterative scheme and a preconditioner using this block structure are proposed. Also, a spectral preconditioner,  which makes use of the information in a Krylov subspace to preconditionate the next system is studied. Second and fourth order exponential methods are also proposed to integrate the time dependent neutron diffusion equation, where the exponential of the system matrix has to be multiplied by a vector. These schemes allow us to work without building explicitly the system matrix, and different methods are compared to calculate the product of the system matrix by a vector, such as a Krylov method, a Chebyshev approximation method, and a Leja Points method.

Some situations arise in which a set of modes of a nuclear reactor core have to be updated, as in perturbative calculations or in the use of modal methods. Updating this set of modes can be very expensive when using the Arnoldi method, and for this reason several methods based on a Newton Iteration are proposed, such as the Modified Block Newton method, a One Sided Block Newton method and a Two Sided Block Newton method. As an alternative strategy, a reduced order model to update a set of modes based on the Proper Generalized Decomposition is also proposed. This method obtains an approximated solution as a sum of separable functions over the whole domain, reducing the multidimensional problem to a set of one-dimensional problems.