How to calculate the energy band diagram for extrinsic and intrinsic semiconductor devices? I am seeking advice for methods and equations to calculate the energy band diagram of intrinsic and extrinsically doped semiconductor devices. For example an InP PIN junction. I am hoping to be able to reproduce a diagram like the one below for a variety of materials at various dopings and structures.

 A: You have four explicit parameters, semiconductor class (intrinsic or doped), intrinsic type (pure or compound), dopant type (for the doped), and dopant level (for the doped). You may have two implicit parameters, temperature and externally applied field.
The staring point is to collect information for the band gaps and dopant levels for the range of materials that you want to consider. These values are published; they are not values that you need to calculate. Many are in standard textbooks for materials science or solid state physics. Harder-core citations can be found in published journals or citation handbooks (e.g. CRC Handbook of Chemistry and Physics).
You now have collected the reference data that you need.
The next step should be to create the pictures at zero external field and O K. Whatever method you choose, by hand or by a computer tool, you will have to "draw on" the table of reference values from above. You have not asked directly for a recommendation for a specific computer package for this part. Perhaps this is a good thing, because recommendations here would likely be personal (subjective), ranging from the "graphical" UI/UX approaches in MatLab, Mathematica, Maple, or Igor Pro to the "roll your own code" approaches in python or C-variant languages.
The final step would be to incorporate the equations that define how band gaps vary with temperature or dopant amounts (e.g. ppm). Here again, searching in textbooks for materials science or solid state physics will reveal the starting points if not the details for the equations that you need. Feed these equations into your hand drawings or code.
If your goal is to produce this as a teaching resource (e.g. to present the dynamics in a classroom demo), I could point to a related example that I have developed to show dynamically how $\ln(\sigma)$ versus $1/T$ graphs change as a function of dopant type and level. Contact me privately or note in the comments that you would appreciate the link.
