NaN in converting a sym to double

I have troubles in converting a sym to double using the functions double, vpa and matlabFunction. I need to solve a system of 6 equations in 6 variables.

Here is the code:

% Here I define the symbolic variables and I assign the physical values to the known constants

clear all
clc
syms s C1 C2 C3 C4 C5 C6

c33=(8.674*10^10); 
rho=(2.65*10^3);
e33=(0.2);
p3=(4*10^-4);
beta3=(1.16*10^6);
CE=(782);
epsilon3=(0.392*10^-10);
k3=(1.4);
R0=(10/rho);
T0=(293);
v=(0.5);
tau=(0.05);
L=(1); 
A=(0.01); 
c0=sqrt(c33/rho);
eta0=rho*CE/k3;

% Here I define other parameters which depends from the known constants

a2=c0*eta0*L*e33^2/(epsilon3*c33);
a3=beta3*T0/c33;
b1=c0*eta0*L;
b2=p3*T0/e33;
d1=beta3/(rho*CE);
d2=p3*c0*e33*L/(k3*epsilon3);
d=rho*R0/v;

%Here i define the parameters which depends on variable s

f0=b2*s^2*(b2*d2-b1)/(a2*b2-a3*b1); 
f2=1-b2*(d1+s/(a3*(1+tau*s)))+a2*b2^2*s/(a3*(a2*b2-a3*b1)*(1+tau*s))-(b2*d2-b1)*(b2+a3+b2*s^2*A^2)/(a2*b2-a3*b1);
f4=(b2/(a3*s*(1+tau*s)))*(1+(s*A)^2-a2*(b2+a3+b2*(A*s)^2)/(a2*b2-a3*b1));
f=d*b2/v

f = 23.4400

lambda1=sqrt((-f2+sqrt(f2^2-4*f4*f0))/(2*f4));
lambda2=-lambda1;%ok
lambda3=sqrt((-f2-sqrt(f2^2-4*f4*f0))/(2*f4));
lambda4=-lambda3;%ok

K=f/(f4*(s/v)^4+f2*(s/v)^2+f0);

j1=(b2*s^2/lambda1^2-a3-b2-(s*A)^2*b2)/(b2*a2-b1*a3);
j2=j1;
j3=(b2*s^2/lambda3^2-a3-b2-(A*s)^2*b2)/(b2*a2-b1*a3);
j4=j3;
j=(b2*v^2-b2-a3-(s*A)^2*b2)/(b2*a2-b1*a3);

l1=lambda1*(b1*j1-1)/b2;
l2=-l1;
l3=lambda3*(b1*j3-1)/b2;
l4=-l3;
l=s*(1-b1*j)/v*b2;

%Here i define the system of 6 equations in the 6 variables C1,C2,C3,C4,C5,C6

eq1=C1+C2+C3+C4+K==0

eq1 = 

eq2=C1*exp(lambda1*c0*eta0)+C2*exp(-lambda1*c0*eta0)+C3*exp(lambda3*c0*eta0)+C4*exp(-lambda3*c0*eta0)+K*exp(-s*c0*eta0/v)==0

eq2 = 

eq3=j1*(C1+C2)+j3*(C3+C4)+j*K+C6==0

eq3 = 

eq4=j1*(C1*exp(lambda1*c0*eta0)+C2*exp(-lambda1*c0*eta0))+j3*(C3*exp(lambda3*c0*eta0)+C4*exp(-lambda3*c0*eta0))+j*K*exp(-s*c0*eta0/v)+C5*c0*eta0+C6==0

eq4 = 

eq5=l1*lambda1*(C1+C2)+l3*lambda3*(C3+C4)-l*K*(s/v)==0

eq5 = 

eq6=l1*lambda1*(C1*exp(lambda1*c0*eta0)+C2*exp(-lambda1*c0*eta0))+l3*lambda3*(C3*exp(lambda3*c0*eta0)+C4*exp(-lambda3*c0*eta0))-l*K*(s/v)*exp(-s*c0*eta0/v)==0

eq6 = 

%Here I solve the system using solve and assign the solutions to the
%variables C1Sol, C2Sol, C3Sol etc.

sol=solve([eq1,eq2,eq3,eq4,eq5,eq6],[C1,C2,C3,C4,C5,C6])

sol = struct with fields:
    C1: (293*(57778692707576729461947844999582017804524826614641173017171257660101514430426230006279123831546424940719604492187500*exp(-(19995904568807233200750958508269*((5646...
    C2: -(293*(57778692707576729461947844999582017804524826614641173017171257660101514430426230006279123831546424940719604492187500*exp((19995904568807233200750958508269*((5646...
    C3: -(293*(57778692707576729461947844999582017804524826614641173017171257660101514430426230006279123831546424940719604492187500*exp(-(19995904568807233200750958508269*((564...
    C4: (293*(57778692707576729461947844999582017804524826614641173017171257660101514430426230006279123831546424940719604492187500*exp((19995904568807233200750958508269*((56469...
    C5: -(28123333861010285370486549721925529650597577007769742311686144*(266307650218891705634070777582211948957102890095849377951258158030818029652256063063962731685726415988...
    C6: (11910695183454101875630379991929377555328*(1790345147005850324723858317095972251906830842615772176737675478990582619998576332772812925860663697832281675536809253536572...

C1Sol=sol.C1;
C2Sol=sol.C2;
C3Sol=sol.C3;
C4Sol=sol.C4;
C5Sol=sol.C5;
C6Sol=sol.C6;

The CiSol (i=1,...,6) solutions i get from this system are symbolic expressions dependent on s. When i try to evaluate some of them using subs and then double:

D=subs(C2Sol,{s},{1})

D = 

G=double(D)

G = NaN

what i get as output is NaN. I get an actual result only if I use number <10^(-2), but i need this functions to work with s ranging from 5 to 30.This happens also if i use vpa or matlabFunction instead of double.

I know that this procedure (or something like this) was used in matlab, since it is reported in some scientific papers (doi:10.1088/0964-1726/18/2/025003 and doi: 10.1016/j.heliyon.2018.e00860). They solve the same system and then use the solutions to define functions (u, phi, theta) which they use to perfom numerical Laplace anti-trasform. I cannot do this, since the CiSol that i get seems to not work very well.

Thank you in advance for any help or suggestion.

Answer 1

[Torsten]

exp(lambda1*c0*eta0)

and

exp(lambda3*c0*eta0)

evaluate to Infinity for s big enough (s slightly bigger than 0 already suffices).

[Raffaele]

Thanks a lot, I actually have doubts about the c0*eta0 factor.