This article considers the optimal control of nuclear fission reactors modeled by parabolic partial differential equations. The neutrons ape divided into fast and thermal groups with two equations describing their interaction and fission, while a third equation describes the temperature in the reactor. The coefficient for the fission and absorption of the thermal neutron is assumed to be controlled by a function through the use of control rods in the reactor. The object is to maintain a target neutron flux shape, while a desired power level and adjustment costs are taken into consideration. A nonlinear optimality system of six equations is deduced, characterizing the optimal control. An iterative procedure is shown to contract toward the solution of the optimality system in small time intervals. The theory is extended to include the effect of other fission products, leading to coupled ordinary and partial differential equations. Numerical experiments ape also included, suggesting directions for further research.