I think you have to use numerial method of integration such as simpson rule or trapezoidal rule. Make a routine for either of these methods and calculate them for each z values and add up. If sum over them, you will get the answer. I hope it works for you. Let me know if you have problem on making subroutine.
Goof luck.
