Schrodinger equation numerical solution

[aaaa202]

I have written a program that solves the Schrödinger equation numerically using the finite difference method based on the attached article. The end goal is to make a program that solves the Schrödinger and Poisson equations self-consistently for the conduction band in different heterostructures.
To do so you need to invoke different relations. For example that the occupation of the kth eigenenergy is given by:

n_k = m*/(πħ2) ∫E_k∞ 1/(1+exp((E-E_f)/kT)) dE

, where E_k is the fermi energy and m* is the effective mass in the band.
Now the problem I have come to is that my units create problems in the above equation. Using the electron mass for the effective mass (I don't know what I should use, I guess it depends on the material) and ħ and a trial potential in the same scale as these units I get eigenenergies of size ≈10-34.
But since kT≈10^(-25) for the systems I am working with the exponential in the above equation simply yields 1 for all energies.
I feel like an amateur that this is giving me problems, but I really don't know what to do at this point. Should I use a different effective mass? I guess I should but is that the main problem? Should I add some kind of conduction band offset? And also what fermi energy should I use for my system?
Hope you will take time to answer some of my questions :)

 

[mfb]

The electron mass should not appear in the exponential. E-E_f is of the order of electronvolts.

 

[aaaa202]

But surely E_k should depend on the electron mass?

 

[aaaa202]

The trouble is that the values of E close to E_k inserted in the exponential simply yields 1. But maybe that isn't a problem actually.
Another thing is however: How do I know which value of the fermi energy to use?

 

[mfb]

I'm not sure how you defined E_k - the energy levels depend on the electron mass but there are no extremely small or large values appearing so the calculations should not run into problems.
Also, you seem to mix E_k and E_f in your posts here.

The trouble is that the values of E close to E_k inserted in the exponential simply yields 1.

So what? Energy levels very close to the Fermi energy will be filled by 1/2, that is the correct result.

How do I know which value of the fermi energy to use?

If the Fermi energy is wrong, the total number of electrons will be too low or too high. You can iterate (guess energy -> calculate number of electrons -> improve guess -> ...) until you find the right value, but there might be some trick to directly estimate it if the energy levels allow that.