matlab fatlab?

So I am at work again and I am going to write a few thoughts: Phase can be calculated from reflectivity using the the integral

phi(E)= -1/pi  int 0 to infinity   ln R(X/E)  / (X^2 -1) dX   (this is understood to be a principal part integral)

So the main issues for me are as follows:

1.  We only know the function R for a limited range so we need to extrapolate it for other values of energy (frequency).  For the moment, I am just assuming that R is constant outside that range. I do not know if this a reasonable assumption or how to go about saying what the reflectivity is for large frequencies.

For large frequency, I get the error as  Log(R_infinity)/(2 Pi) Log |  ( Emax-E)/(Emax+E) |
For small frequency, Log(R_0)/(2 Pi)   Log |  (E-Emin)/(E+Emin)|

R_infinity = R(E_max)
R_0=R(E_min)

2. I have to find an appropriate method of discretizing the integral that takes into account that we need to take a principal part integral.
2a. For the moment, I am using a right-side riemann sum. Perhaps I will implement simpson's rule. The disadvantage of the latter method is that it assumes equal spacing so I cannot see if the large spacing is the issue.
2b. Say we want to estimate the phase for a frequency omega. Then we need to know the reflectivity very precisely near omega. However, from what I hear from the lab, the reflectivity is measured every 0.2 eV, which I do not think is accurate enough.

Anyways, so I spent most of today trying to get the phase to predict more accurately. I had an idea that would let me integrate over the singularity. Basically, I just wanted to subtract the singularity out around a small symmetric neighborhood. This contributes nothing to the principal part integral but removes the singularity. This seemed fine but when I went to try and integrate this new function, I spent probably 3 hours trying to figure out how to get matlab to integrate this.  This was pretty frustrating for me because I could have done this in 15 mins in mathematica but my lab does not use that program. I understand that this is just something that I will learn in time but I don't like having work go so slowly like this. I eventually asked one of the guys that I work with a few questions that I was able to figure it out. The main idea is that many functions need to be able to handle a vector input. I also didn't realize that when you write equations, you need to allow for vectors to work as well.