Looks like you are just getting the exact solutions in rational form a/b
, the evaluation to floating point should match the numerical integrations
fn_sym = double(fn)
fn_num = areaOfRect - quadl(@(p)(1 - (1 - p.^5).^10),0.01,threshold)
fp_sym = double(fp)
fp_num = quadl(@(p)(1 - (1 - p.^5).^10),threshold,1.01)