The noninteracting electronic structures of tight binding models on bipartite lattices with unequal numbers of sites in the two sublattices have a number of unique features, including the presence of spatially localized eigenstates and flat bands. When a emph{uniform} on-site Hubbard interaction $U$ is turned on, Lieb proved rigorously that at half filling ($rho=1$) the ground state has a non-zero spin. In this paper we consider a `CuO$_2$ lattice (also known as `Lieb lattice, or as a decorated square lattice), in which `$d$-orbitals occupy the vertices of the squares, while `$p$-orbitals lie halfway between two $d$-orbitals. We use exact Determinant Quantum Monte Carlo (DQMC) simulations to quantify the nature of magnetic order through the behavior of correlation functions and sublattice magnetizations in the different orbitals as a function of $U$ and temperature. We study both the homogeneous (H) case, $U_d= U_p$, originally considered by Lieb, and the inhomogeneous (IH) case, $U_d eq U_p$. For the H case at half filling, we found that the global magnetization rises sharply at weak coupling, and then stabilizes towards the strong-coupling (Heisenberg) value, as a result of the interplay between the ferromagnetism of like sites and the antiferromagnetism between unlike sites; we verified that the system is an insulator for all $U$. For the IH system at half filling, we argue that the case $U_p eq U_d$ falls under Liebs theorem, provided they are positive definite, so we used DQMC to probe the cases $U_p=0,U_d=U$ and $U_p=U, U_d=0$. We found that the different environments of $d$ and $p$ sites lead to a ferromagnetic insulator when $U_d=0$; by contrast, $U_p=0$ leads to to a metal without any magnetic ordering. In addition, we have also established that at density $rho=1/3$, strong antiferromagnetic correlations set in, caused by the presence of one fermion on each $d$ site.