The short answer is, yes, computational chemistry gets them wrong because we don't have a way to calculate orbital energies exactly. In the end, the root cause is the general problem with quantum many-fermion systems - we can't solve the Schrödinger equation exactly.
Let's have a look to different computational chemistry methods to understand why that's the case. In general, more sophisticated methods will reduce the energy gap between alpha and beta electrons, but as we will see, none will guarantee that it doesn't exist.
Hartree-Fock methods
Hartree-Fock methods are inherently approximate: when we solve a Hartree-Fock calculation, we are solving not a system of N electrons, but N coupled systems of electrons that only interact indirectly. Of the energy terms in the Hamiltonian, HF accounts for the kinetic energy of each electron, its interaction with all the nuclei, the coulombic repulsion between electrons, and the exchange energy with other electrons with the same spin (because they are fermions), but we ignore electronic correlation. In HF, electrons only interact with each other indirectly - by contributing to a background external potential in the resolution of each electron's particular Hamiltonian. This makes it necessary to solve the equations iteratively using a mean field approach - what we usually call "self-consistent field" or SCF. Correlation requires explicitly considering more than one electron at the same time, so it's not compatible with the core HF model.
But, even if it's approximate, shouldn't we still get the same, not exact, values for alpha and beta electrons? - sure we do, in systems where the error introduced by the HF assumptions is the same for alpha and beta electrons.
For instance, if you run an unrestricted HF calculation on singlet $\ce{O_2}$, where you have the same number of alpha and beta electrons, you'll get exactly the same values for the alpha and beta orbitals (and these will be exactly the same as you'd get for a restricted HF calculation of the same system).
However, for triplet $\ce{O_2}$, alpha and beta electrons are exposed to slightly different interactions - since correlation depends on the spin of the interacting electrons, ignoring it affects alpha and beta electrons differently.
Post-Hartree-Fock methods
Post-Hartree-Fock methods attempt to correct the HF limitations by first resolving a Hartree-Fock system and then adding a correction for electronic correlation. The problem, however is that the corrections aren't exact.
Møller-Plesset perturbation methods take the HF wavefunction and then add a numerical correction to the energy in a perturbative fashion - you apply the HF method to double-excited (MP2), triple-excited (MP3), and so on versions of the same ground state HF wavefunction and then introduce a correction to the energy of your orbitals (not to the wavefunctions themselves) based on that. It has several problems - you are only accounting for dynamic correlation (how electrons influence the movement of other electrons within your HF wavefunction), and, because excited estates are not variationally limited, there's no guarantee that you will not overestimate your correction and will end up with an energy below the exact energy. But even if those problems could be completely corrected, since part of the correlation is not accounted for, it will still give slightly different energies for alpha and beta orbitals in a system with $\mathrm{S > 0}$.
Configuration interaction (CI) methods apply the HF method not only to the ground state but also to a number of excited states, and then find the correct solution to the system variationally, as a linear combination of those configurations. This method is in principle exact - and it's variational, so we are guaranteed not to underestimate the energy - the problem is that it is a series expansion, in which we add a term for each possible configuration state function (CSF). Still, even a full CI taking into account all the possible CSFs won't account for all of the correlation - CI is only exact in the limit of a complete basis set. In practice, expanding the series is computationally very expensive, and it is truncated relatively early - full CI is only practical for the simplest of systems.
When it's applicable, however, difference in energy between alpha and beta electrons should be negligible.
Coupled cluster
Coupled cluster (CC) has a theoretical foundation that is distinct from HF and post-HF methods; essentially, it takes a wavefunction (usually a Slater determinant) and performs a number of excitations on it, producing a coupled cluster description of the system that, like full CI, is in principle exact. However, like CI, CC is a series expansion - you have terms for single, double, triple, etc. excitations, and we would require all possible excitations to be taken into account to fully consider electron correlation. As with CI, we usually truncate the series quite early - the usual "gold standard" of computational chemistry, against which all other methods have usually been benchmarked in the last couple of decades, is CCSD(T) - meaning that we consider the first and second term, and the third term only partially. This is usually considered to be as exact as it gets in computational chemistry, but again, it is prohibitive for all but the simplest of systems.
Like with full CI, CC differences between alpha and beta electrons will be negligible.
Density functional theory
I left DFT for the end because it's an entirely different beast, theoretically speaking, but in the end the problem comes down to the same issue: you don't account for correlation correctly, and that flaw impacts alpha and beta electrons differently for $\mathrm{S > 0}$.
Density functional theory, in principle, shouldn't have the same problems Hartree-Fock-related methods have, since we only with a single, global, three-dimensional density function $\mathrm{\rho(r)}$; and, as per the second Hohenberg-Kohn theorem, the ground state electron density is unique and variationally limited, so in principle we could keep applying the functional to different three-dimensional functions until we find the absolute minimum - and that's our exact ground state density.
There, are, however, two big problems.
First, in practical terms, we need a way to translate the electron density to wavefunctions so we can understand what's happening at the electron level. For instance, the very values for orbital energies that you mention for a DFT calculation are not inherently in the density function that DFT deals with; there's no way to distinguish the first from the second occupied orbitals in the density function. What we do is introduce a formalism to be able to translate the density function to orbital terms - that's the Kohn-Sham formalism, which boils down to treating the global particle density function, which distributes all N particles in space, as the sum of N single-particle functions, each corresponding to one electron, and which interact with each other only indirectly through a background external potential, which accounts for the coulomb, exchange and correlation interactions, so we update densities and potentials iteratively in another mean field, self-consistent approach. You'll notice that this is pretty much the same thing we did with the wavefunction in Hartree-Fock.
Two considerations, however. First, the Kohn-Sham formalism allows us to map the total density to orbital occupations and from there calculate a number of properties that we couldn't derive from the total density function, but it guarantees only that the sum of the single-electron densities is equal to the total density function for the ground state - not that the depiction of a particular orbital is accurate. So, values for specific molecular orbitals (such as orbital energy) should be taken with a pinch of salt. In fact, in practice, orbital properties from DFT are most accurate around the HOMO-LUMO gap, while values for core orbitals or virtual orbitals can diverge quite wildly from reference values. This is so because orbitals are constructed in a way that best fits density, and not one that best fits the wavefunction. You can find a lot of discussions about whether Kohn-Sham orbitals are physically meaningful and how they relate to the actual orbitals of the system.
As you can imagine, if you have important spin effects in your system, spin deviation in these properties are almost guaranteed.
But there's also an advantage over Hartree-Fock: in HF, we had to drop correlation because we couldn't incorporate it - it required treating more than one electron wavefunction explicitly. But, in DFT, wavefunctions are an afterthought - what we are minimising is the total density function, so there's nothing that prevents us from introducing correlation. We could even ditch the Kohn-Sham formalism and calculate correlation directly in terms of the total density function. So, aside from the problem of breaking down the total density function into a number of one-electron densities, we could in principle at least get an exact total alpha electron density and a total beta electron density, couldn't we?
Unfortunately, no, because we don't know how the functional for electrons looks like. And, if we move from wavefunctions to density, we also lose the ability to calculate exchange interactions exactly. The only system for which we can calculate these exactly is the free electron gas; so all the functionals we regularly use are proposals for the exchange or correlation terms, or both, that approximate the exchange-correlation interactions but that don't represent them exactly.
And as with all the other cases, if we are not capturing all the exchange and correlation effects, that will impact alpha and beta electrons differently in a system with $\mathrm{S > 0}$.
