How can I make this model that iterates over 1000 times without converging more efficient?

Question

I have this code that previously had simpler versions of p and q functions that worked. I have modified them for my model. So I know that the problem arised when I changed these. I have checked so many times and cannot figure out if I made a mistake and I have almost given up hope.

Parameters:

(* Cost share of housing capital in housing production function *)
\[Beta] = 0.6;
(* Scaling on housing production function *)
g = 0.0005;
(* Radians available for construction *)
(* benchmark \[Theta] is 3 *)
\[Theta] = 2*Pi;
(*elasticity of income*)
\[Gamma] = 1;
(*price elasticity of housing*)
\[Epsilon] = 1;
(*constant in demand function*)
\[Omega] = 1;

Functions:

p[x_, y_, t_, ta_, f_, fa_, 
   u_] := ((\[Epsilon] - 
       1) ((\[Gamma] - 
           1) u + ((y - (f + 
              fa) - (t - ta) x))^(1 - \[Gamma]))/((\[Gamma] - 
          1) \[Omega]))^(-1/(\[Epsilon] - 1));
q[x_, y_, t_, ta_, f_, fa_, 
   u_] := \[Omega] ((\[Epsilon] - 
        1) ((\[Gamma] - 
            1) u + ((y - (f + 
               fa) - (t - ta) x))^(1 - \[Gamma]))/((\[Gamma] - 
           1) \[Omega]))^(\[Epsilon]/(\[Epsilon] - 
        1)) ((y - (f + fa) - (t - ta) x))^\[Gamma];
S[x_, y_, t_, ta_, f_, fa_, 
   u_] := (1/(p[x, y, t, ta, f, fa, u] (\[Beta]) (g)))^(1/(\[Beta] - 
       1));
r[x_, y_, t_, ta_, f_, fa_, 
   u_] := (p[x, y, t, ta, f, fa, 
       u] (g)) ((1/(p[x, y, t, ta, f, fa, 
           u] (\[Beta]) (g)))^(\[Beta]/(\[Beta] - 1))) - 
   1 ((1/(p[x, y, t, ta, f, fa, u] (\[Beta]) (g)))^(1/(\[Beta] - 1)));
h[x_, y_, t_, ta_, f_, fa_, 
   u_] := (g) S[x, y, t, ta, f, fa, u]^(\[Beta]);
Density[x_, y_, t_, ta_, f_, fa_, u_] := 
  h[x, y, t, ta, f, fa, u]/q[x, y, t, ta, f, fa, u];
L[xbar_, y_, t_, ta_, f_, fa_, u_] := \[Theta]*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*\ 
     Density[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
pbar[xbar_, y_, t_, ta_, f_, fa_, u_, 
   pop_] := (\[Theta]/pop)*(1/30000)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*
     Density[x, y, t, ta, f, fa, u]*
     p[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
qbar[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*
     Density[x, y, t, ta, f, fa, u]*
     q[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
hbar[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*
     Density[x, y, t, ta, f, fa, u]*
     h[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
xhat[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x^2*
     Density[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
rev[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := 
  pop*(fa + ta*xhat[xbar, y, t, ta, f, fa, u, pop]);

Case-specific parameters:

Population = 800000;
y = 70000;
t = 350;
f = 3000;
ta = 250;
fa = 1000;
ra = 45000;

Running the model:

Lu[u_] := (solution = 
    FindRoot[{r[xbar, y, t, ta, f, fa, u] - ra == 0}, {xbar, 10}, 
     AccuracyGoal -> 1]; xbar /. solution);

fu[u_] := N[L[Evaluate[Lu[u]], y, t, ta, f, fa, u] - Population];

findRoot[fun_, lower_, upper_, eps_] := 
  Module[{ul = lower, uh = upper, count = 0}, 
   While[And[Abs[fun[ul]] > eps, count < 1000], 
    If[fun[(uh + ul)/2] > 0, ul = (uh + ul)/2, uh = (uh + ul)/2]; 
    count = count + 1]; Print["Iterations: ", count,
    " Population: ", Population, " people",
    " utility: ", N[ul],
    " City boundary (xbar): ", Lu[N[ul]], " km",
    " Average house price (pbar): ", 
    pbar[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " euros",
    " Price at 5 km from the centre p(5): ", 
    p[5, y, t, ta, f, fa, N[ul]]/30000, " euros",
    " Price at 10 km from the centre p(10): ", 
    p[10, y, t, ta, f, fa, N[ul]]/30000, " euros",
    " Average density: ", 2*Population/(\[Theta]*Lu[N[ul]]^2), 
    " people per ?",
    " Average floor to area ratio hbar: ", 
    hbar[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population],
    " Density at 5 km from the centre Density(5): ", 
    Density[5, y, t, ta, f, fa, N[ul]], " people per ?",
    " Density at 20 km from the centre Density(20): ", 
    Density[20, y, t, ta, f, fa, N[ul]], " people per ?",
    " Average distance to the centre xhat: ", 
    xhat[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " km",
    " Tax revenue from transport: ", 
    rev[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " euros", 
    Plot[{Density[x, y, t, ta, f, fa, N[ul]]}, {x, 0, Lu[N[ul]]}]]]; 
findRoot[fu, 3000, 10000, 1];

After running this (which took a whole night), all values were at an extreme, that didn't make sense. I am not sure how that happened, but I was hoping I could get some help figuring out.

I think the answer lies in p and q functions if that helps anything. It worked fine before I changed those, but is it just too much for my computer to handle and should I increase the iterations or is there something I'm not seeing. Any help would be greatly appreciated.