Convert the symbolic functions to function handles, and use numerical integrals -
P_l=50;
v=0.1;
k=15;
Tm=1;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
nu=0.3;
syms x z x_prime z_prime t dT_dx dT_dz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-((sqrt((x-v*t)^2 + (z)^2) + (x-v*t)))/(2*alpha)))/(4*3.14*k*sqrt((x-v*t)^2 + (z)^2)) +T0;
dT_dx=diff(T,x);
dT_dx_prime=subs(dT_dx,[x,z],[x_prime,z_prime]);
dT_dz=diff(T,z);
dT_dz_prime=subs(dT_dz,[x,z],[x_prime,z_prime]);
Gxh= (1/(4*pi))*(3*(xm/(xm^2 + zp^2)) + 2*(xm*zm^2/(xm^2 +zm^2)^2))-(1/pi)*(3*(xm*(z_prime*zp + xm^2)/(xm^2 + zp^2)^2)-(3*(z_prime)^2*xm*zp*2 +xm^3*(4*z_prime^2 + 6*z*z_prime + z^2 + xm^2))/(xm^2+zp^2)^3);
Gxv= (-1/(4*pi))*((zp/(xm^2 + zp^2))+ 2*((xm^2*zm/(xm^2+zm^2)^2)-(xm^2*zm)/(xm^2 +zm^2)^2))-(1/(2*pi)*(2*(zp/(xm^2 + zp^2)))-((2*z-z_prime)*(zp^2-xm^2)/(xm^2+zp^2)^2)+(2*z*z_prime*zp*(3*xm^2-zp^2))/(xm^2+zp^2)^3);
p= P_l*exp(-(-v*t)/2*alpha)/(4*pi*k*(-v*t));
%Convert to a function handle
T = matlabFunction(T);
p = matlabFunction(p);
fun = matlabFunction(Gxh * dT_dx + Gxv * dT_dz);
%Define terms as funciton handles
term1 = @(t,x,z) integral2(@(x_prime,z_prime) fun(t,x,z,x_prime,z_prime),-inf,inf,0,inf);
term2 = @(x,z) integral(@(t,x,z) (2.*z)./pi*(p(t).*(t - x).^2./((t - x).^2 + z^2).^2), t, -inf, inf);
term3 = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*v);
Sigma = term1(0,0.001,0) + 0 + term3(0,0.001,0)
