MATLAB: How to solve for limit of double integral

double integralfzerolimit

I am new to MATLAB. Please help me. I want to solve for the limit of double integral.
syms u x
a = 5;
b = 0.2;
c = 10;
d = 0.1;
k=0.01;
m = 1./(2.^(a+c-1).*gamma(a).*gamma(c).*b.^a.*d.^c);
h(u,x) = m.*(u+x).^(a-1).*(u-x).^(c-1).*exp(-(u+x)./(2.*b)-(u-x)./(2.*d));
f(x) = int(h,u,abs(x), inf);
solPk = fzero(@(Pk) integral(f,-inf,Pk)-k,[-10,10]);
Thank you so much.

Best Answer

solPk = vpasolve(int(f,x,-inf,Pk) == k)
There is a closed form integral, but MATLAB is not able to find it. It is
piecewise(Pk < 0, -(770558495002939/5159780352000000000000)*exp(10*Pk)*(75937500*Pk^9-296156250*Pk^8+601425000*Pk^7-817897500*Pk^6+808258500*Pk^5-594641250*Pk^4+322528500*Pk^3-123369750*Pk^2+29996190*Pk-3512131), ...
0 <= Pk, 1310719999999999239/1310720000000000000+(1/25798901760000000000000)*(-11650844444444437680000*Pk^4-40389594074074050624000*Pk^3-58409566814814780902400*Pk^2-41590925590123432642560*Pk-12267389871934149256192)*exp(-5*Pk))
You could fzero() or fsolve() on that after making a guess about whether Pk will be positive or negative.