The problem of optimally controlling xenon spatial oscillations is formulated as a problem in the calculus of variations for distributed parameter systems. The resulting partial differential equations (space- and time-dependent) are then approximated by a nodal representation to obtain a set of ordinary differential equations (time-dependent) with mixed (initial and final) boundary conditions. An iterative solution scheme, which utilizes a quasilinearization of the equations and a transformation matrix relating initial to final values of certain variables, is employed to obtain numerical results. Feasibility of the method is established by several sample calculations. A physical interpretation is given the Lagrange multiplier functions which initially are introduced for mathematical considerations.