Grant, I "think" I managed to vectorize (or more generally "array-cize") the functions. There is a weird tweak on the CND function that solved the problem that put me in a quandary for a bit. Using pylab library (matplotlib) I made the price an array and the 2 strikes of the bull call spread an array for this go around. The trick here is that both arrays needed to be the SAME dimension where it is len(price) by len(strikes) or in this example 2 rows by 16 columns. For the price ranges are repeated on each row as such: [[20.,21,.....,30.,31.,32.,33.,34.,35.], [20.,21,.....,30.,31.,32.,33.,34.,35.]] and for the strikes they are as follows (remember same dimensions as the price array): [[25.,25.,......,25.,25.,25], [30.,30.,.......,30.,30.,30.]] If you are wondering why the decimal point, that just makes it float, it does not work if they are integers in the formulas. Since my head is starting to hurt, I decided not the vectorize the time and volatility variables at this stage, but I think it generalizes if you keep track of the indices. The array just gets bigger and keeps adding more dimensions. I felt for the the different expirations, I just looped it 3 times was easier. Hindsight being 20/20, I think looping the function for all practical purposes is fast enough, unless you are looking to attempt to use this realtime in a quick fuse daytrading tool....NOT. It was just a good exercise for me to go through to understand the matplotlib/pylab interface with python. I attached (hope it works) a png file of a bull call spread for your review.