function [p_n, rho_n] = solve_general_balance_equations(lambda_n, mu_n)
% SOLVE_GENERAL_BALANCE_EQUATIONS - Solve the steady-state balance equations given arrays of \lambda_n and \mu_n 
% 
% Input: 
% 
% lambda_n = [ \lambda_0, \lambda_1, \lambda_2, \cdots ] 
% mu_n     = [ \mu_1,     \mu_2,     \mu_3,     \cdots ]  
% 
% Written by:
% -- 
% John L. Weatherwax                2007-09-10
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

rho_tmp = lambda_n ./ mu_n;           % <- do pointwise division ... 
rho_n   = [ 1, cumprod( rho_tmp ) ]; 

% compute the normailizing constant: 
p_0 = 1/sum( rho_n );

% compute the steady state probabilites: 
p_n = p_0*rho_n;