spainMath.max

lp: fullratsimp(
	(N-x-y)/B - h*y/(1-h*B) - x/(1-B)
)$

n: fullratsimp(num(lp))$

d: denom(lp)$

write("Denominator is positive (≥ a positive number): ")$
fullratsimp(subst(h=1-ε, d - B*(B-1)^2));

write("Numerator is an upward-pointing quadratic, so first solution is the maximum (negative slope of parabola is the second derivative of likelihood)")$
nb: collectterms(expand(n), B);
sol: solve(nb=0, B)$

write("Numerator is still positive at B=0, so first solution is positive")$
fullratsimp(subst([B=0], n));

write("Numerator is negative at B=1, so second solution is out of bounds (≥1)")$
fullratsimp(subst([B=1], n));

string(sol[1]);