As already pointed out in a comment by John Custer, Dirac was 100% correct in
The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation.
For light systems, typically thought as H-Kr, non-relativistic quantum mechanics i.e. the Schrödinger equation is sufficient to describe chemistry; for heavier nuclei you need to use relativistic quantum mechanics i.e. the Dirac equation which is a bit more complicated. In many cases we can invoke the Born-Oppenheimer method, and assume that the nuclei are moving in the instantaneous field of the electrons; now all we need to solve is the electronic problem.
We know that the exact solution to the electronic problem in this case is achievable with the full configuration interaction (FCI) a.k.a. exact diagonalization method, in which you describe the electronic many-particle wave function by a weighted sum of electronic configurations a.k.a. Slater determinants as $|\Psi \rangle = \sum_k c_k |\Phi_k\rangle$. These electronic configurations are built by distributing $N$ electrons into $K$ single-particle states a.k.a. orbitals. For the method to be accurate, $K\gg N$, and actually to get the exact solution you need $K \to \infty$.
Now, to find the ground state (as well as any excited states), you just need to diagonalize the many-electron Hamiltonian in the basis of the electron configurations. But, the problem is that the number of the electron configurations grows extremely rapidly.
If we assume that we are looking at a spin singlet state, then you have $N/2$ spin-up and $N/2$ spin-down electrons. For each spin, there are ${K \choose N/2} = \frac {K!} {\frac N 2 ! (K-\frac N 2)!} $ ways to populate the orbitals. The total number of electron configurations for $N$ electrons in $K$ orbitals, typically denoted as ($N$e,$K$o), is then ${K \choose N/2}^2 = \left[ \frac {K!} {\frac N 2 ! (K-\frac N 2)!} \right]^2$.
Even for the case of a very small number of orbitals, $K=N$, the number of configurations quickly becomes huge. (8e,8o) has 4900 configurations, (10e,10o) has 63 504, (12e,12o) has 853 776, (14e,14o) has 11 778 624, (16e,16o) has 165 636 900, (18e,18o) has 2 363 904 400, (20e,20o) has 34 134 779 536, and (22e,22o) has 497 634 306 624.
Although you can still solve the (8e,8o) problem with dense matrix algebra on modern computers, you see that very quickly you have to become very smart in how to diagonalize the matrix. Because the Hamiltonian is a two-particle operator, it is extremely sparse in the basis of the electron configurations: if two configurations differ by the occupation of more than two orbitals, the Hamiltonian matrix element is zero by Slater and Condon's rules. Moreover, for large problem sizes you also want to avoid storing the matrix, which is why you want to use an iterative method. (The famous Davidson method for iterative diagonalization was actually developed exactly for this purpose!)
With smart algorithms, billion-configuration calculations i.e. the (18e,18o) problem have been possible since the early 1990s, see e.g. Chem. Phys. Lett. 169, 463 (1990). However, despite the huge increase in computational power in the last 30 years, the barrier has barely budged: as far as I am aware, the largest FCI problem that has been solved is the (22e,22o) calculation in J. Chem. Phys. 147, 184111 (2017).
The thing to note here is that even the (22e,22o) calculation is not large enough to solve a single atom in an exact fashion: you need a lot more orbitals to achieve quantitative accuracy with experiment. Although a high-lying orbital contributes only very little to the correlation energy, there are A LOT of them.
Exactly as Dirac wrote, approximations are needed. Density-functional approximations are extremely popular in applications, but they are far from exact. On the other hand, high-accuracy studies often employ the coupled-cluster method, which is a reparametrization of the FCI method; however, it, too would exhibit exponential scaling were it not truncated - i.e. approximated.