We consider two steady state heat conduction systems called, S and Sα, in a multidimensional bounded domain D for the Poisson equation with source energy g. In one system we impose mixed boundary conditions (temperature b on the boundary Γ1, heat flux q on Γ2 and an adiabatic condition on Γ3). In the other system, the condition on Γ1 is replaced by a convective heat flux condition with coefficient α. For each of these systems, we consider three associated optimization problems (Pi) and (Piα), i = 1, 2, 3 where the variable will be the source energy g, the heat flux q and the environmental temperature b, respectively. In the particular case that D is a rectangle, the explicit continuous optimization variables and the corresponding state of the systems are known. In the present work, by using a finite difference scheme, we obtain the discrete systems (Sh) and (Sh α) and discrete optimization problems (Pih) and (Pihα), i = 1, 2, 3, where h is the space step in the discretization. Explicit discrete solutions are found, and convergence and estimation errors results are proved when h goes to zero and when α goes to infinite. Moreover, some numerical simulations are provided in order to test theoretical results. Finally, we note that the use of a three–point finite–difference approximation for the Neumann or Robin boundary condition at the boundary improves the global order of convergence from O(h) to O(h2).