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

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;
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;
eq3=j1*(C1+C2)+j3*(C3+C4)+j*K+C6==0;
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;
eq5=l1*lambda1*(C1+C2)+l3*lambda3*(C3+C4)-l*K*(s/v)==0;
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;

%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]);
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})
G=double(D)
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.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

% 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

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

%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]);
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})
G=double(D)
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. 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

% 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

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

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

Cart

Cart