Fortran program called pmodesz and some of its subroutines which was used to solve the mode equation for the modes that are well behaved at the horizon. Note that the program calls a subroutine odebs which comes from the book Numerical Recipes and which therefore cannot be given here. But any differential equation solver could be used to replace it.
Here are more subroutines: convert, metric, pseries0, pseries

Fortran program called qmodesz and some of its subroutines which was used to solve the mode equation for the modes that are well behaved at infinity. Note that the program calls a subroutine odebs which comes from the book Numerical Recipes and which therefore cannot be given here. But any differential equation solver could be used to replace it.
It calls the convert and metric subroutines above and also qseries

Output from the pmodesz and qmodesz programs is read into the program pqwkb8z. This program uses the convert subroutine above. It also uses the subroutine wkb8. It computes the products of the p and q modes and their derivatives that are necessary for the stress-energy tensor. It also normalizes those products correctly. It then subtracts off the 8th order wkb approximation for these products.

Output from pqwkb8z is read into the program lsumz where the 5 sums over the angular momentum parameter ell are computed for each value of the frequency omega. No subrouines are called by lsumz.

Output from lsumz is read into the final program omintz where the integral over the frequency omega is computed for each of the
five sums over ell. Since the renormalization process involves subtracting the 8th order wkb approximation in pqwkb8z, it must be added back on and the
renormalization counterterms are subtracted from it. A high frequency expansion is also added and subtracted. See Phys. Rev. D51, 4337 (1995) for details
The subroutine convert above is called. Also called are the subroutines Tmna and Traceanom.

The data files for the curve in Figure 1 with A33 = B33 = 0 are: one, two, and three.
In each file the first column is the value of the radial coordinate s = (r-r0)/r0 with r0 the radius of the event horizon. The second column is , the third is Ttt, the fourth is Trr. For the stress-energy tensor components one index is up and one is down. For the rest refer to the file omintz above or send an email message to Paul Anderson.
For A33 = 0, B33 = 1 the data files are: one, two, and three.
In the first one a quantity related to one of the mode sums is in the third column, Ttt is in the fourth column, and Trr is in the fifth column. In the other two and all of the rest of the data files Ttt is in the third column and Trr is in the fourth column.
For A33 = 0, B33 = 2 the data files are: one, two, and three.
For A33 = B33 = 2, the data files are: one and two.

The data files for Figure 2 are the ones listed above for A33 = B33 = 0 and for A33 = 0, B33 = 2.

The data files for the curves in Figure 3 with A33 = B33 = 0 are one and two.
The data files for the curves with A33 = 4, B33 = 6 are one and two
The columns are the same as those for most of the data files for Figure 1.