I think the issue is that your function is not continuous at $H= \frac{Y-\delta Y}{\alpha +r}$ as at that value of $H$ there is a division by zero.
If $H\leq \frac{Y-\delta Y}{\alpha +r}$, then g1[H]=1 and h1[H]=0. Otherwise, g1[H]=1 and h1[H]=0. That gives two functions:
f1 = ((1 + (1 - p) β) ((-H r + Y)^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/
(1 - 1/σ) + (p β ((Y - H (r + α))^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/
(1 - 1/σ)
f2 = (p β ((Y δ)^(1 - 1/ϵ) + 5000^(1 - 1/ϵ) θ^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ) +
((1 + (1 - p) β) ((-H r + Y)^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ)
Maximize both functions (with the appropriate restrictions on $H$), and then take the maximum result:
Y = 1000;
β = 0.95;
θ = 10;
α = 0.05;
r = 0.05;
ϵ = 1.2;
σ = 0.2;
δ = 0.3;
p = 0.51;
sol1 = Maximize[{f1, H <= (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0}, H]
(* {-3.26868*10^-24, {H -> 7000.}} *)
sol2 = Maximize[{f2, H > (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0}, H]
(* {-4.45483*10^-24, {H -> 14852.5}} *)
max = If[sol1[[1]] > sol2[[1]], sol1, sol2]
(* {-3.26868*10^-24, {H -> 7000.}} *)
This approach seems more stable when additional restrictions are added:
sol1 = Maximize[{f1, H <= (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0, H < 14000}, H]
(* {-3.26868*10^-24, {H -> 7000.}} *)
sol2 = Maximize[{f2, H > (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0, H < 14000}, H]
(* {-4.46379*10^-24, {H -> 14000.}} *)
max = If[sol1[[1]] > sol2[[1]], sol1, sol2]
(* {-3.26868*10^-24, {H -> 7000.}} *)