Getting imaginary values as solutions to equations using solve function

I am trying to solve these three equations in matlab. I am getting the answer as a function of z and imaginary values after substituting the values of z and a warning. Where am I going wrong..can I get only real valued solutions? Also, I want to plot sol1 as a function of z. Can someone please help me with this?

P_l=50;
v=0.1;
k=15;
Tm=1375;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
h= 100*10^9;
G= 150*10^9;
nu=0.3;
psi= 1- exp(-(0.45*h)/(2*G));
syms x z x_prime z_prime t dT_dx dT_dz sigma_xx sigma_yy sigma_zz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-v*((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]);
%Define Green's function
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));
Gzh= (-1/(4*pi))*(3*(xm/(xm^2 +zp^2)- (xm/(xm^2 +zm^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));
Gxzh= (1/(4*pi))*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))+2*((zp*xm^2/(xm^2 +zp^2)^2)-(zm*xm^2/(xm^2 +zm^2)^2)))-(1/pi)*(3*(z_prime*zp^2 +xm^2*(2*z+z_prime)/2*(xm^2+zp^2)^2)-(z_prime^3*(z_prime^2 + 3*z*z_prime + 3*z^2) + z^3*z_prime^2 + xm^2*(z_prime^3 + 6*z*z_prime^2 + 6*z^2*z_prime + z^3)+z*xm^4)/(xm^2 +zp^2)^3);
Gzv= (1/(4*pi))*(3*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))) +2*((zm*xm^2/(xm^2 +zm^2)^2)-(zp*xm^2/(xm^2 + zp^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));
Gxzv= (xm/(4*pi))*((1/(xm^2+zp^2) - 1/(xm^2+zm^2)) + 2*((zp^2/(xm^2 +zp^2))-(zm^2/(xm^2+zm^2)^2)))-(xm/(2*pi))*((4*z*zp/(xm^2 +zp^2)^2) + (2*z*z_prime*zp*(3*zp^2 -xm^2)/(xm^2 + zp^2)^3));
%Convert to a function handle
T = matlabFunction(T);
p = matlabFunction(p);
fun1 = matlabFunction(Gxh * dT_dx + Gxv * dT_dz);
fun2 = matlabFunction(Gzh * dT_dx + Gzv * dT_dz);
fun3 = matlabFunction(Gxzh * dT_dx + Gxzv * dT_dz);
%Define terms as funciton handles for Sigma_xx_therm
term1_xx = @(t,x,z) integral2(@(x_prime,z_prime) fun1(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xx = @(x, z) (2 * z) / pi * integral(@(t) (p(t) .* (t - x).^2 ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xx = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xx_therm = -(alpha*E/(1-2*nu))*term1_xx(0,0,-0.0005) + term2_xx(0,-0.0001) + term3_xx(0,0,-0.0005);
%Define terms as funciton handles for Sigma_zz_therm
term1_zz = @(t,x,z) integral2(@(x_prime,z_prime) fun2(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_zz = @(x, z) (2 * z^3) / pi * integral(@(t) (p(t) ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_zz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_zz_therm = -(alpha*E/(1-2*nu))*term1_zz(0,0,-0.0005) + term2_zz(0,-0.0001) + term3_zz(0,0,-0.0005);
%Define terms as funciton handles for Sigma_xz_therm
term1_xz = @(t,x,z) integral2(@(x_prime,z_prime) fun3(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xz = @(x, z) (2 * z^2) / pi * integral(@(t) (p(t) .* (t - x)./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xz_therm = -(alpha*E/(1-2*nu))*term1_xz(0,0,-0.0005) + term2_xz(0,-0.0001) + term3_xz(0,0,-0.0005);
% components of deviatoric stress
S_xx=(2*sigma_xx - sigma_yy - sigma_zz)/3;
S_yy= (2*sigma_yy - sigma_xx - sigma_zz)/3;
S_zz=(2*sigma_zz - sigma_yy - sigma_xx)/3;
S_eq=sqrt(sigma_xx^2 + sigma_yy^2 + sigma_zz^2);
n_xx=S_xx/S_eq;
n_yy=S_yy/S_eq;
n_zz=S_zz/S_eq;
%Equations
eqn1= (1/E)*(sigma_xx-nu*(sigma_yy+sigma_zz))+ (1/h)*(sigma_xx*n_xx + sigma_yy*n_yy + sigma_zz*n_zz)*n_xx == psi*((1/E)*(Sigma_xx_therm - nu*(sigma_yy + Sigma_zz_therm))+(1/h)*(Sigma_xx_therm*n_xx + sigma_yy*n_yy + Sigma_zz_therm*n_zz)*n_xx);
eqn2= (1/E)*(sigma_yy-nu*(sigma_xx+Sigma_zz_therm))+(1/h)*(sigma_xx*n_xx + sigma_yy*n_yy +sigma_zz*n_zz)*n_yy==0;
eqn3= sigma_yy==0.5*(sigma_xx+sigma_zz);
sol=solve([eqn1,eqn2,eqn3],[sigma_xx,sigma_yy,sigma_zz]);
sol1=sol.sigma_xx;
sol2=sol.sigma_yy;
sol3=sol.sigma_zz;

[Dyuman Joshi]

If you want to use solve, it would be better to stick to using symbolic expressions/functions, and avoid converting to function handles.

Answer 1

[John D'Errico]

P_l=50;
v=0.1;
k=15;
Tm=1375;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
h= 100*10^9;
G= 150*10^9;
nu=0.3;
psi= 1- exp(-(0.45*h)/(2*G));
syms x z x_prime z_prime t dT_dx dT_dz sigma_xx sigma_yy sigma_zz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-v*((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]);
%Define Green's function
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));
Gzh= (-1/(4*pi))*(3*(xm/(xm^2 +zp^2)- (xm/(xm^2 +zm^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));
Gxzh= (1/(4*pi))*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))+2*((zp*xm^2/(xm^2 +zp^2)^2)-(zm*xm^2/(xm^2 +zm^2)^2)))-(1/pi)*(3*(z_prime*zp^2 +xm^2*(2*z+z_prime)/2*(xm^2+zp^2)^2)-(z_prime^3*(z_prime^2 + 3*z*z_prime + 3*z^2) + z^3*z_prime^2 + xm^2*(z_prime^3 + 6*z*z_prime^2 + 6*z^2*z_prime + z^3)+z*xm^4)/(xm^2 +zp^2)^3);
Gzv= (1/(4*pi))*(3*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))) +2*((zm*xm^2/(xm^2 +zm^2)^2)-(zp*xm^2/(xm^2 + zp^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));
Gxzv= (xm/(4*pi))*((1/(xm^2+zp^2) - 1/(xm^2+zm^2)) + 2*((zp^2/(xm^2 +zp^2))-(zm^2/(xm^2+zm^2)^2)))-(xm/(2*pi))*((4*z*zp/(xm^2 +zp^2)^2) + (2*z*z_prime*zp*(3*zp^2 -xm^2)/(xm^2 + zp^2)^3));
%Convert to a function handle
T = matlabFunction(T);
p = matlabFunction(p);
fun1 = matlabFunction(Gxh * dT_dx + Gxv * dT_dz);
fun2 = matlabFunction(Gzh * dT_dx + Gzv * dT_dz);
fun3 = matlabFunction(Gxzh * dT_dx + Gxzv * dT_dz);
%Define terms as funciton handles for Sigma_xx_therm
term1_xx = @(t,x,z) integral2(@(x_prime,z_prime) fun1(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xx = @(x, z) (2 * z) / pi * integral(@(t) (p(t) .* (t - x).^2 ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xx = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xx_therm = -(alpha*E/(1-2*nu))*term1_xx(0,0,-0.0005) + term2_xx(0,-0.0001) + term3_xx(0,0,-0.0005);
%Define terms as funciton handles for Sigma_zz_therm
term1_zz = @(t,x,z) integral2(@(x_prime,z_prime) fun2(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_zz = @(x, z) (2 * z^3) / pi * integral(@(t) (p(t) ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_zz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_zz_therm = -(alpha*E/(1-2*nu))*term1_zz(0,0,-0.0005) + term2_zz(0,-0.0001) + term3_zz(0,0,-0.0005);

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 9.8e+16. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

%Define terms as funciton handles for Sigma_xz_therm
term1_xz = @(t,x,z) integral2(@(x_prime,z_prime) fun3(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xz = @(x, z) (2 * z^2) / pi * integral(@(t) (p(t) .* (t - x)./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xz_therm = -(alpha*E/(1-2*nu))*term1_xz(0,0,-0.0005) + term2_xz(0,-0.0001) + term3_xz(0,0,-0.0005);
% components of deviatoric stress
S_xx=(2*sigma_xx - sigma_yy - sigma_zz)/3;
S_yy= (2*sigma_yy - sigma_xx - sigma_zz)/3;
S_zz=(2*sigma_zz - sigma_yy - sigma_xx)/3;
S_eq=sqrt(sigma_xx^2 + sigma_yy^2 + sigma_zz^2);
n_xx=S_xx/S_eq;
n_yy=S_yy/S_eq;
n_zz=S_zz/S_eq;
%Equations
eqn1= (1/E)*(sigma_xx-nu*(sigma_yy+sigma_zz))+ (1/h)*(sigma_xx*n_xx + sigma_yy*n_yy + sigma_zz*n_zz)*n_xx == psi*((1/E)*(Sigma_xx_therm - nu*(sigma_yy + Sigma_zz_therm))+(1/h)*(Sigma_xx_therm*n_xx + sigma_yy*n_yy + Sigma_zz_therm*n_zz)*n_xx);
eqn2= (1/E)*(sigma_yy-nu*(sigma_xx+Sigma_zz_therm))+(1/h)*(sigma_xx*n_xx + sigma_yy*n_yy +sigma_zz*n_zz)*n_yy==0;
eqn3= sigma_yy==0.5*(sigma_xx+sigma_zz);
sol=solve([eqn1,eqn2,eqn3],[sigma_xx,sigma_yy,sigma_zz]);
sol.sigma_xx

ans = 

vpa(sol.sigma_xx)

ans = 

So it looks like three roots, Two of which are complex conjugates. Your problem is probably implicitly a cubic polynomial, so that should not be a surprise.

vpa(sol.sigma_yy)

ans = 

vpa(sol.sigma_zz)

ans = 

If you want only the real roots, then keep only root number 1. Where is the problem?

As far as plotting a solution as a function of z, since z is not present in the solution, that does not make complete sense.