function r = eq_return_v09(rho,Parameter)
   
    Solution0    = zeros(3,1);
    Solution0(1) = log((Parameter.alpha/(1-Parameter.alpha)) ...
                   /(Parameter.delta+rho+Parameter.sigma*Parameter.g)); 
             % log of normalized capital (k in notes)
    Solution0(2) = Parameter.govbonds;                                          
             % normalized bond holdings of households (b in notes)
    Solution0(3) = log(rho+(Parameter.sigma-1)*Parameter.g);            
             % log of return to bonds (rb in notes) - g
             
    
    
    % - load solver -----------------------------------------------------------
    options=optimset('disp','off','LargeScale','off', ...
                 'TolFun',.00000001,'MaxFunEvals',10^7,'MaxIter',10^7);
    Solution = fsolve(@(Solution)equations_v09(Solution,Parameter),Solution0,options);
    Check=equations_v09(Solution,Parameter);
    assert(max(abs(Check))<0.001); 
    
    k_norm  = exp(Solution(1));
    rb      = Parameter.g+exp(Solution(3));
    r       = (Parameter.alpha/(1-Parameter.alpha))*(rb-Parameter.g)/k_norm - Parameter.delta;
  
end
