Question) How should I solve this equation from a finite-state Markov chain?

I have 2,001 states, and say the transition matrix of Markov chain is given. Then there is a variable x(s) attached to each state s \in \{1,...,2001\}, and I have to solve the following equation for any x(s):

\frac{a}{f(x(s))}+x(s)=b+\mathbb{E}_{s}[\frac{1}{f(x(s'))}] .

where the function f is non-linear.

I understood this problem as solving a system of non-linear equations, so I used nlsolver package.
I could solve for 50 states, but for 2,001 states, I think I will never get my result as I have to simulate this problem 10,000 times.

This is a replication exercise, and the author mentions that this is a fairly simple numerical exercise. Would there be another approach to solve this problem more time efficiently?

So you are looking for a fixed point of x=T(x').
If T is a contraction mapping, iteration over x=T(x') will converge to the fixed point.
But this approach requires that you can efficiently compute x(s) from the LHS of your equation. That depends on what f looks like.

Do you have any additional information about x (e.g. monotonicity)?

Thank you very much!

I could solve with NLsolver package fixedpoint function.
It takes only like 1/10 of time now!