自己贴出答案，希望能帮到大家！
clear all; close all
F = @(x) 1/2*erf(x/sqrt(2));
xs = [0:0.05:1].'; % 变量x的采样序列
ys = [10000:5000:100000].'; % 变量y的采样序列
[XS, YS] = meshgrid(xs, ys);
ZS = 0*XS;
for g = 1:length(xs)
    x = xs(g);
    for  k = 1:length(ys)
        y = ys(k);
        xy = x*y;
        ZS(k, g) = 1.104*xy + 0.72*quad(F,0,xy/40);
    end
end
mesh(XS,YS,ZS);

呵呵。自娱自乐型的。

不是，自己慢慢弄出来了！