Any good book in Semiconductor Physics will have a description of the k.p method. Try Fundamentals of Semiconductor Physics by Peter Yu and Manuel Cardona. Another reference for Kane Model and EFA are chapters 2 and 3 of "Wave Mechanics Applied to Semiconductor Heterostructures" by Gerald Bastard.
If you want a more mathematically/group theory oriented discussion, try the first chapters on the extended Kane model in "Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems" by Roland Winkler.
But the rationale is very simple and it's a simple consequence of Bloch's theorem, that states that the eigenstates of an electron sitting in a periodic potential can be written as:
$\psi_{k}(\vec{x}) = e^{i\vec{k}\cdot\vec{x}} u_{\vec{k}}(\vec{x})$
where, this is very important, $u(\vec{x})$ have the same periodicity of the potential.
Now, you can just plug this identity in Schrödinger's equation:
$\left(\frac{p^2}{2m} + V(\vec{x})\right)\psi_{\vec{k}}(\vec{x}) = E(\vec{k}) \psi_{\vec{k}}(\vec{x})$
and you will eventually find a differential equation for the $u_{\vec{k}}(\vec{x})$ functions(*):
$\left(\frac{p^2}{2m} + \frac{\hbar}{m} \vec{k}\cdot\vec{p} + V(\vec{x})\right)u_{\vec{k}}(\vec{x}) = E(\vec{k}) u_{\vec{k}}(\vec{x})$
Notice that this looks like an eigenvalue problem for a $\vec{k}$-dependent hamiltonian. For each value of $\vec{k}$, you can solve this problem and find the eigen-energies. This will give rise to the energy bands. Notice the term $\vec{k}\cdot\vec{p}$ in this hamiltonian, that gives the method its name. It looks like this k-dependent hamiltonian is just the old hamiltonian plus a perturbation term proportional to $\vec{k}\cdot\vec{p}$. (**)
Now, you can construct a matrix representation of this eigenvalue problem by projecting it on a basis. For a convenient choice of basis you must specify the symmetry of the specific material you are dealing with. You are going to study the irreducible representations of the symmetry group of the crystal.
The trick is that the function $u_{\vec{k}}(\vec{x})$ must be invariant under the same symmetry group as the original hamiltonian
$H = \frac{p^2}{2m} + V(\vec{x})$
If this is the case, then you can easily choose a basis in which this hamiltoniean is diagonal! You just have to choose vectors that transform according to the irreducible representations of this group. Now you only have to worry about the k.p term. The matrix hamiltonian wil be a diagonal part plus an off-diagonal perturbation. But this off-diagonal term is not arbitrary. If you study how the operator $\vec{k}\cdot\vec{p}$ transforms relatively to the various irreps of the group, you'll have lots of zeros. You can determine all of the terms in the hamiltonian that will be non-zero with this procedure.
Now you just have a matrix with a bunch of constants in each entry where symmetry tells you that the hamiltonian shouldn't in principle be zero. Now you have to find this constants. You have two possibilities: you get data from experiments, or from ab initio simulations. This is as far as you can go by just using symmetry.
Note that until this point this is all exact. No approximations were made. The next step from the Kane model is to include the possibility of slightly broken translational symmetry. This takes us to the Envelope Function Approximation (EFA). But this is already getting too long and confuse! :)
Please take a look in the Yu-Cardona book. There's a really nice, pedagogical discussion about this there. Bastard's book have a good derivation of the Kane matrix and the EFA matrix. Wrinkler's book have a more mathematically rigorous derivation from group theory arguments.
If you need experimentally measured constants for your Kane matrix there's a very comprehensive review by Vurgaftman et al for III-V semiconductors:
Band parameters for III-V compound semiconductors and their alloys. J. Appl. Phys., v. 89, n. 11, p. 5815–5875, 2001.
(*) if you have spin-orbit interactions or other relativistic effects this will complicate a bit the derivation but it will remain essentially valid.
(**) There's a crucial difference: this new eigenvalue problem have different boundary conditions! You can solve it just inside a unit cell, cause $u(\vec{x})$ is periodic. And the boundary condition is given by this periodicity requirement.
