No, you’re missing a ./R
on the second term. This should be it:
using OrdinaryDiffEq, Plots
#Parameters
P = 75.0 #pression 50 => collapse
#initial conditions
R₀ = [5.0]
dR₀ = [1.0]
tspan = (0.0,50.0)
#numerical
function rp(ddR, dR, R, P, t)
@. ddR = (P - (1.5 *dR^2))/R
end
#pass to the solver
probleme = SecondOrderODEProblem(rp, 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)
savefig("RP1.pdf")
Fixing that term seems to stabilize the solution a lot more too, so I don’t think a method for stiff equations is necessary.