MATLAB ODE15s: imaginary parts

I have to solve $$\left\{\begin{matrix} \ddot{R}=-\frac{3}{2}\frac{\dot{R}^{2}}{R}+\frac{1}{R\rho}\left [P_{v}(T_{\infty})-P_{0}+\frac{(L\rho_{g})^{2}}{hT_{\infty}}\dot{R}-4\mu \frac{\dot{R}}{R}-\frac{2S}{R} \right ]\\ h=\frac{\lambda}{2R}\left [ 2+\left ( \frac{-6\rho c_{p}\dot{R}}{\pi h}\right )^{1/3}-\left ( \frac{12\rho c_{p}\dot{R}}{\pi h} \right ) \right ]\\ \ln\left ( \frac{\rho_{g}}{\rho_{c}}\right )=c_{1}\theta^{1/3}+c_{2}\theta^{2/3}+c_{3}\theta^{4/3}+c_{4}\theta^{3}+c_{5}\theta^{37/6}+c_{6}\theta^{71/6} \end{matrix}\right.$$ where $\theta=1-\frac{1}{T_{c}}\left ( T_{\infty}+\frac{L\rho_{g}}{h}\dot{R} \right )$ and $R(0)=R_{0}$, $\dot{R}(0)=0$. If a use a pressure $P_{0}=2100$, the code works; if I set $P_{0}=2000$, MATLAB returns this error: "Warning: Imaginary parts of complex X and/or Y arguments ignored."

This is my code:

function out = densita_20_latente(t,y,R0)
S = 0.073; %N/m
rho = 998;
mi = 1.005e-3;
P0 = 2100; %<-----------------------------------------
L = 2454000;
Tinf = 293;
lambda = 0.5982;
cp = 4184;
Pv_Tinf = 2329.6;
c1 = -2.03150240;
c2 = -2.6830294;
c3 = -5.38626492;
c4 = -17.2991605;
c5 = -44.7586581;
c6 = -63.9201063;
Tc = 647.096;
rhoc = 322;
%
out=[y(2)
-1.5 * y(2)^2 / y(1) + 1 / rho / y(1) * (Pv_Tinf - P0 +...
(L * y(4))^2 / Tinf * y(2) / y(3) -...
4 * mi * y(2) / y(1) - 2 * S / y(1))
-y(3) + lambda * 0.5 / y(1) * (2 + (-6 * rho * cp * y(2) / pi / y(3))^(1/3)...
- 12 * rho * cp * y(2) / pi / y(3))
-y(4) + rhoc * exp(c1*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(1/3)...
+ c2*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(2/3)...
+ c3*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(4/3)...
+ c4*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(3)...
+ c5*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(37/6)...
+ c6*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(71/6)) ];


run file

R0=510e-6;
y0=[510e-6; 0; 1172.94; 0.01716];
M = [1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0];
options = odeset('Mass',M);
[t,y] = ode15s(@(t,y) densita_20_latente(t,y,R0),[0 46.75e-6],y0,options);
t = t*1000000;
y = y*1000000;
[t,y(:,1)];
plot(t,y(:,1))


Edit: I am going to simplify the problem I have to solve:

$$\left\{\begin{matrix} \ddot{R}=-\frac{3}{2}\frac{\dot{R}^{2}}{R}+\frac{1}{R\rho}\left [P_{v}(T_{\infty})-P_{0}+\frac{(L\rho_{g})^{2}}{hT_{\infty}}\dot{R}-4\mu \frac{\dot{R}}{R}-\frac{2S}{R} \right ]\\ h=\frac{\lambda}{2R}\left [ 2+\left ( \frac{-6\rho c_{p}\dot{R}}{\pi h}\right )^{1/3}-\left ( \frac{12\rho c_{p}\dot{R}}{\pi h} \right ) \right ]\\ \end{matrix}\right.$$

function out = convezione_20(t,y,R0)
S = 0.073;
rho = 998;
mi = 1.005e-3;
P0 = 2000; %<-------------------------------------------
L = 2454000;
Tinf = 293;
rho_g = 1/57.79;

lambda = 0.5982;
cp = 4184;
Pv_Tinf = 2329.6;
%
out=[y(2)
-1.5 * y(2)^2 / y(1) + 1 / rho / y(1) * (Pv_Tinf - P0 +...
(L * rho_g)^2 / Tinf * y(2) / y(3) -...
4 * mi * y(2) / y(1) - 2 * S / y(1))
- y(3) + lambda * 0.5 / y(1) * (2 + (-6 * rho * cp * y(2) / pi / y(3))^(1/3)...
- 12 * rho * cp * y(2) / pi / y(3)) ];


run file

R0 = 510e-6;
y0 = [R0; 0; 1172.94];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M);
[t,y] = ode15s(@(t,y) convezione_20(t,y,R0),[0 40e-6],y0,options);
t = t*1000000;
y = y*1000000;
[t,y(:,1)];
plot(t,y(:,1))


Same problem as before. With 2100 the code works; with 2000 MATLAB reads "Warning: Imaginary parts of complex X and/or Y arguments ignored." I tried using as exponent $1$ instead of $1/3$ in the expression of $h$, but the warning now becomes "Warning: Failure at t=2.874099e-007. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.021085e-021) at time t."

• The error message that you're getting indicates that your function is returning complex values rather than real values. Most likely this is because some of the quantities being raised to fractional powers are negative. Nov 22 '16 at 16:04
• @BrianBorchers I edited the original post. According to what you said: can the problem be with $(-6\rho c_{p}\dot{R}/(\pi h))^{1/3}$? Nov 22 '16 at 17:04
• Are there some time derivatives missing on the left side, for instance in the equation for $h$? Is the equation for $\rho_g$ "algebraic" or (when corrected) differential? Nov 27 '16 at 14:56