If there any possible solution (analytical/numerical) of a this nonlinear third order differential equation? The equation I'm trying to solve is the following
$$m_3\ddot{y}^3 + m_1\ddot{y} + k_1 y + k_3y^3 = f_1(t) $$
All the parameters are real $m_1,m_3,k_1,k_3\in \mathbb{R}$ and constant.
$f_1(t)\in \mathbb{R}$ is a function of $t$ as is $y(t)$.
Obviously if $m_3=0$ then the solution can be found because the equation takes the form
$$\ddot{y} =  a y + by^3 + c,$$
where $a=-\dfrac{k_1}{m_1}$, $b = -\dfrac{k_3}{m_1}$, $c = \dfrac{f_1}{m_1}$,
which is the classical form of a nonlinear ODE as
$$ \ddot{y} = F(y,t).$$
EDIT:
Please find below my implementation in MATLAB of DAE with NO success so far.
First I define the DAE function as
function out = funTest(t,y,yp)

% coefficients
m1 = 9.8845614144e-10;
m3 = 2.00030511601069e-12;

f1 = 1.66666666666667e-05;

k1 = 0.214188819381172;
k3 = 8.6610733212416;

omega = 2*pi*(-82.454*t + 10030);

out = [yp(1) - y(2)
    yp(2) - y(3)
    m3*y(3)^3 + m1*y(3) + k1*y(1) + k3*y(1)^3 - f1*sin(omega*t)];

Below the main file
tspan = [60 120];
M = [1 0 0; 0 1 0; 0 0 0]; % constant mass matrix
options = odeset('RelTol',1e-4,'Jacobian',{[],M});

y0 = [0; 0; 0];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@funTest,0,y0,[1 1 0],yp0,[],options);

[t,y] = ode15i(@funTest,tspan,y0,yp0,options);


plot(t,y(:,1));
ylabel('y');

I tried changing both options as
options = odeset('RelTol',1e-4,'AbsTol',[1e-10 1e-6 1e-6], ...
   'Jacobian',{[],M});

but nothing happens
 A: The problem is to avoid the potential complex solutions for the cubic. Using a symbolic processor we can follow with the cubic real root (always exists). Using MATHEMATICA as the symbolic solver we submit
sols = Solve[m3 y''[t]^3 + m1 y''[t] + k1 y[t] + k3 y[t]^3 == f[t],y''[t]]

obtaining as real root 
$$
y''(t) \to \frac{\sqrt[3]{2} \left(9 m_3^2 \left(f(t)-y(t) \left(k_3 y(t)^2+k_1\right)\right)+\sqrt{3} \sqrt{m_3^3 \left(27 m_3
   \left(-f(t)+k_3 y(t)^3+k_1 y(t)\right){}^2+4 m_1^3\right)}\right){}^{2/3}-2 \sqrt[3]{3} m_1 m_3}{6^{2/3} m_3 \sqrt[3]{9
   m_3^2 \left(f(t)-y(t) \left(k_3 y(t)^2+k_1\right)\right)+\sqrt{3} \sqrt{m_3^3 \left(27 m_3 \left(-f(t)+k_3 y(t)^3+k_1
   y(t)\right){}^2+4 m_1^3\right)}}}
$$
and after that we solve numerically the DE
m1 = 9.8845614144 10^-10;
m3 = 2.00030511601069 10^-12;
f1 = 1.66666666666667 10^-05;
k1 = 0.214188819381172;
k3 = 8.6610733212416;
omega = 2*Pi*(-82.454*t + 10030);
f[t_]:= f1*Sin[omega*t]
tmax = 1;
equ = y''[t] == (y''[t]/.sols[[1]]);
soly = NDSolve[{equ, y[0] == 0.1, y'[0] == 0.1},y,{t,0,tmax}][[1]];
Plot[Evaluate[y[t] /. soly], {t, 0, tmax}]


