Understanding, deriving, writing, testing and debugging an ab initio code can be a lenghty and tedious task. I'd like to provide a starting point for you here. If you just to it for pedagocial reasons, it might be advisable to start with the atomic problem and try to solve it with DFT. The effort for that is not too big, but it covers nearly all the nesseary physical principles to understand how one arrives at total energies and electron densities. With such a program, you can predict for example, energy levels of a certain element. Also, If you want to extend your program to treat crystalline material i. e., the written routines for the atomic problem will be of fundamental use. Core states are essential calculated like this.
The following websites provide some beginner-friendly introduction to the topic of simple DFT calculations of atoms.
I want to list down some the steps that are nessary for a simple solver written in python.
Setting up the equations
The starting point is the Kohn–Sham equation, which is the one-electron Schrödinger-like equation of a fictitious system of non-interacting particles. If you take the Born-Oppenheimer approximation into account and assume that the Kohn-Sham potential is radial symmetric, the task reduces to solve the radial part of the Kohn-Sham equation under some restrictions. Basically, you need to construct a solver that calculates the radial wave function and the eigenenergy for a given Kohn-Sham potential. For this you will need a shooting method similar to . Basically, one solves the an inital value problem starting from $r=0$ and $r=\infty$ and change the energy till the solution matches at some point.
Calculating the Kohn-Sham potential
The effective Kohn-Sham potential is the sum of the core contributions ($1/r$ coloumb potential), the exchange-correlation (input here is the electron density) and the Hartree potential. In order to get the Hartree potential, you will need to solve the radial Poisson equation.
In order to set up a self-consistent cycle, you will need to calculate the poisson equation, which can be treaten as an inital value problem or also as an boundary value problem.
Remark to computational efficency in Python
The computational effort is considerable for bigger atoms. You might want to take a look some acceleration methods for Python like Numba.