How can I solve for total money spent if I use a non-constant MPC?
I don't know that the distribution you claim is wrong but it strikes me as having a couple of very undesirable features. For one, virtually the entire population all of the has >0.99 MPC, which seems very high for a distribution you think is centered around 3/4. For another, if you think that the income shocks in your model are transitory rather than permanent, my reading of the literature is that these are way too large. But that's all separate from your question about estimating things conditional on the actual MPC distribution and ordering. 

Ordering from low to high MPC seems to take a strong stance on how trade works in the economy. Again, I don't know that is wrong but it doesn't strike me as true. If we think that MPC is caused by income then this is akin to saying the richest guy get the money first, spends some it at the second richest guy's store and saves the rest, he spends the spending of guy 1 at the third richest guy and so on. The economy strikes me as less neat than that and so a random ordering makes more sense to me in the absence of evidence about where high vs. low MPC people spend their money. 

I tried messing around with simulations (in Python) to see how well we could approximate these random chains of consumption. In general, experimenting a few parameter values with uniform, beta, and your distribution above, distribution doesn't matter so much in determining total consumption effect. Using the geometric or arithmetic means of the MPC distribution give nice approximations of actual total consumption effect from the random chains of consumption.* On the other hand, if you order the households by MPC from high to low or low to high this makes a huge difference in the resulting total consumption effects. 

* - That said, I explicitly ruled out distributions of MPC where MPC can be zero (or less). If that happens the geometric mean will be zero and the resulting approximation using the geometric mean quite poor. 

Code: 

# Python 3.x code for stack exchange question
# http://economics.stackexchange.com/questions/3123/how-can-i-solve-for-total-money-spent-if-i-use-a-non-constant-mpc
import numpy as np

# Parameters
households = 5001
sims = 10000
initial_shock = 5

# Alternative marginal propensity to consume populations
MPC_vec_lin = np.linspace(0.01, 0.2, households)
MPC_vec_beta = np.random.beta(2, 5, size=households)

# THis is the function Gragas asks for
nvec = np.arange(0,households)
MPC_vec_gragas = (nvec + 2) / np.sqrt(nvec**2 + 5*nvec + 7)

MPC_vec = MPC_vec_gragas

# Calculate geometeric and arithmetic mean MPC from the population
approx_geo_mean = np.mean(np.log(MPC_vec))
approx_long_run_eqiv_MPC = np.exp(approx_geo_mean)
avg_MPC = np.mean(MPC_vec)

# Assuming that all households had the same MPC (@ geo / arithmetic level), what is the total effect?
approx_total_shock_with_geo = initial_shock / (1 - approx_long_run_eqiv_MPC)
approx_total_shock_with_mean = initial_shock / (1 - avg_MPC)

print('Initial shock size:', initial_shock)
print('Approximate total effect with geometric mean:', approx_total_shock_with_geo)
print('Approximate total effect with arithmetic mean:', approx_total_shock_with_mean)

# Simulating many possible orderings, what is the total effect for each ordering?
total_shock_vec = np.zeros(sims)
for i, item in enumerate(total_shock_vec):
    MPC_vec_random_order = np.random.choice(MPC_vec, replace=False, size=households)
    random_products = np.cumproduct(MPC_vec_random_order)
    total_shock = np.sum(initial_shock * random_products) + initial_shock
    total_shock_vec[i] = total_shock
# Compare resulting average and std. of total effect from simulations
print('Average total effect from simulations:',  np.mean(total_shock_vec))
print('Std Dev of total effect from simulations:',  np.std(total_shock_vec))

# Conclusion, mean and geomean do a very nice job of approximating the resulting total effect without knowing the
# general specific ordering. However, the worst and best total effect ordering (low to high and high to low) are
# quite different.

MPC_vec_low_to_high = MPC_vec.copy()
MPC_vec_low_to_high.sort()

MPC_vec_high_to_low = MPC_vec_low_to_high.copy()
MPC_vec_high_to_low = MPC_vec_high_to_low[::-1]

product_low_to_high = np.cumproduct(MPC_vec_low_to_high)
total_shock_low_to_high = np.sum(initial_shock * product_low_to_high) + initial_shock

product_high_to_low = np.cumproduct(MPC_vec_high_to_low)
total_shock_high_to_low = np.sum(initial_shock * product_high_to_low) + initial_shock

print('Maximum effect ordering effect size:', total_shock_high_to_low)
print('Minimum effect ordering effect size:', total_shock_low_to_high)


Results: 

Initial shock size: 5
Approximate total effect with geometric mean: 6023.35708426
Approximate total effect with arithmetic mean: 6154.37897874
Average total effect from simulations: 5999.77721767
Std Dev of total effect from simulations: 642.715838738
Maximum effect ordering effect size: 16675.8496051
Minimum effect ordering effect size: 777.157515589