Matlab double integration result does not match with my self calculate

Question

Here is a double integration, self learning $$ P=\int_{-w}^{w}\int_{l}^{\frac{y_h(x_b+w)}{x_h}+l}\frac{1}{2}\operatorname{erfc}\left[\frac{\log{\frac{z_h(y_b-l)}{y_h}}-\mu}{\sigma\sqrt2}\right]\space dy_bdx_b$$

after derive the final formula. I used matlab to verify my result by codes

fun = @(x,y) (0.5*erfc((log(zh*(y-l)/yh)-mu)/(sqrt(2)*sigma)));
ymax = @(x) (yh/xh*(x+w))+l;
q = integral2(fun,-w,w,l,ymax);

but the result is not compatible with my own integration result. and I think is there something wrong with the MATLAB CODES?