Trying to solve a second order differential equation. Your advices

Hello dear all,

Continuing my work on the Rayleigh-Plesset equation, I wrote this code:

using OrdinaryDiffEq, Plots

#SI units
#external parameters
Pₕ = 1.01325e5                  #hydrostatic pressure (Pa)

#Parameters (liquid = water)
Pᵥ = 2.5e3                      #vapour pressure (Pa)
σ = 73.0e-3                     #surface tension (N/m)
ρ = 1000.0                      #volumic mass (Kg.m-3)
η = 1.0e-3                      #dynamic viscosity (Pa.s)

#Parameter (gas = Ar)
γ = 1.666                       #polytropic ratio (-)

#Acoustic Parameters
Pₐ = 1.3e5                      #acoustic pressure (Pa)
fₐ = 2.0e4                      #acoustic frequency (s-1)

#Initial conditions
R₀ = [4.0e-6]                   #initial radius (m)
dR₀ = [1.0e-2]                  #initial radial velocity (ms-1)
tspan = (0.0,150.0e-6)          #(s) -- 1 acoustic cycle = 50 μs


#Num
function rpnnp(ddR, dR, R, p, t)
    @. ddR = (1/ρ*R)*((Pₕ-Pᵥ-(2*σ/R₀))*(R₀/R)^(3*γ)-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end

P = t->Pₐ*sin(2*pi*fₐ*t + pi/2)   #external acoustic wave

probleme = SecondOrderODEProblem(rpnnp, dR₀, R₀, tspan, P)
solution = solve(probleme, DPRKN6())

@show solution.t
@show solution.u

#graph
plot(solution, vars=(2), linewidth=2,
        xaxis = "t ", yaxis = "R",
        framestyle = :box)

I am expecting a result like in the 4th graph concerning the “radius versus the time” curve:

Rayleigh-Plesset_R-t

However with the code as above the bubble radius evolves hardly.
Do you have a suggestion since I do not see where my error is.
Thank you in advance,
Thierry