Not like this is hard.
Make a loop over k. Generate the roots for each k. Then sort them, taking the middle one in each case.
Or if that is too hard...
syms x k
P = k-3*k*x^2+3*x^3;
r = solve(P,x);
r1 = matlabFunction(vpa(r(1)))
r1 =
@(k)k.*3.333333333333333e-1+k.^2.*1.0./(k.*-1.666666666666667e-1+sqrt((k.*1.666666666666667e-1-k.^3.*3.703703703703704e-2).^2-k.^6.*1.371742112482853e-3)+k.^3.*3.703703703703704e-2).^(1.0./3.0).*1.111111111111111e-1+(k.*-1.666666666666667e-1+sqrt((k.*1.666666666666667e-1-k.^3.*3.703703703703704e-2).^2-k.^6.*1.371742112482853e-3)+k.^3.*3.703703703703704e-2).^(1.0./3.0)
r2 = matlabFunction(vpa(r(2)))
r2 =
@(k)k.*3.333333333333333e-1+k.^2.*1.0./(k.*-1.666666666666667e-1+sqrt((k.*1.666666666666667e-1-k.^3.*3.703703703703704e-2).^2-k.^6.*1.371742112482853e-3)+k.^3.*3.703703703703704e-2).^(1.0./3.0).*(-5.555555555555556e-2+9.622504486493763e-2i)+(k.*-1.666666666666667e-1+sqrt((k.*1.666666666666667e-1-k.^3.*3.703703703703704e-2).^2-k.^6.*1.371742112482853e-3)+k.^3.*3.703703703703704e-2).^(1.0./3.0).*(-5.0e-1-8.660254037844386e-1i)
r3 = matlabFunction(vpa(r(3)))
r3 =
@(k)k.*3.333333333333333e-1+k.^2.*1.0./(k.*-1.666666666666667e-1+sqrt((k.*1.666666666666667e-1-k.^3.*3.703703703703704e-2).^2-k.^6.*1.371742112482853e-3)+k.^3.*3.703703703703704e-2).^(1.0./3.0).*(-5.555555555555556e-2-9.622504486493763e-2i)+(k.*-1.666666666666667e-1+sqrt((k.*1.666666666666667e-1-k.^3.*3.703703703703704e-2).^2-k.^6.*1.371742112482853e-3)+k.^3.*3.703703703703704e-2).^(1.0./3.0).*(-5.0e-1+8.660254037844386e-1i)
The problems is, these roots cross over at smaller k (there is a double root for SOME value of k) so if you ALWAYS want the middle root, then you need to choose the proper branch. (High school math to locate the value of k that results in a double root.)
Finally, you COULD convert the polynomial to a 3x3 matrix, then compute the eigenvalues in a sequence as a function of k. Then use my eigenshuffle to do some of the work for you. Not worth it in this case though.
