误差函数导论与erfi函数的matlab编程(erf based on 2 papers and erfi matlab code)

非常需要，谢谢楼主

终于找到关于误差函数的帖子了

syms w1 w2 a1 b1 d p1 p2 E1 E2 e1 e2 D q0 k p0 R r T f a b s h g Cs p w pi m c x y z
a1=0.6,b1=0.8,D=10000,q0=5400,R=0.07,r=0.023,T=60,s=0.5,h=3000,g=2000,Cs=12000,k=0.25,m=0.01,c=0.016,p0=35000
w1=((1-b1)*D+b1*q0)^2/(4*k*p0*q0*(1-b1)*(1+R*T))
w2=(D-q0)/(k*p0*(1+R*T))
d=((1-b1)*D+b1*q0)^2-4*k*w*p0*q0*(1-b1)*(1+R*T)
p1=((1-b1)*D+b1*q0-sqrt(d))/(2*(1-b1)*k)
p2=((1-b1)*D+b1*q0+ sqrt(d))/(2*(1-b1)*k)
f=exp((log(p)-log(p0)-m*T).^2./(-2*c^2*T))./(sqrt(2*pi*T)*c*p)
u=vpa(f,6)
a=w*p0*q0*(R*T-r*T-a1-a1*R*T)
b=s*q0*T+(h-g)*q0-Cs
E1=a1*int(((1-b1)*(-k*p^2+D*p)+b1*q0*p)*u,p,p2,inf)+(a+b)*(int(u,p,p2,inf)+int(u,p,0,w*p0*(1+R*T)))+a1*int(p*q0*u,p,0,w*p0*(1+R*T))+(w*p0*q0*(R-r)*T+b)*int(u,p,w*p0*(1+R*T),p2)
e1=diff(E1,w)
E2=a1*int(((1-b1)*(-k*p^2+D*p)+b1*q0*p)*u,p,(D-q0)/k,p1)+a1*int(((1-b1)*(-k*p^2+D*p)+b1*q0*p)*u,p,p2,inf)+(a+b)*(int(u,p,0,p1)+int(u,p,p2,inf))+a1*int(p*q0*u,p,0,(D-q0)/k)+(w*p0*q0*(R-r)*T+b)*int(u,p,p1,p2)
e2=diff(E2,w)
x=solve(e1,w)
y=solve(e2,w)

我运行matlab里出现了erfi，并且erfi是因为我上面求了e的x^2这种形式的定积分，不是换数据就能去掉的。
其实我就是想对利润函数E1 E2求导得到e1 e2，令导数为0得到里面唯一那个未知数w，但现在matlab对求导后为0这个方程式算不出一个具体的数，总显示empty sym， 它就不识别erfi这个误差的感觉，怎么能求出最后那个数？
我原以为是积分函数f的精度问题，我就加了个u=vpa(f,6)，但还是解不出来，有没有其他软件能解呢？或者怎么让它能解出来，现在老显示结果empty sym。。。

还有我是菜鸟级别的，看见大家在说求解这种定积分用数值方法，我直接把编程放出来了，能不能直接在上面做改动，不然直接给一段程序我看不懂啊，就要最后能求出w这个数来的，而且是实数不要复数那种形式。。。