rf = rfpnk(x,omega0,data)
solv = vpasolve(rf, [-1000*omega0; 1000*omega0])
x = fsolve(@(x)rfpnk(x,omega0,data), soln, options)
xagain = fsolve(@(x)rfpnk(x,omega0,data), [5.18e15 -2.17e16], options)
But initial values 5.1e+15 -2.1e+16 were not good enough for fsolve() to find something it liked. And notice that even though with the 3 digits of precision of the input solution, that fsolve's f(x) is not great compared to what it is when using the full precision found by solve(): if you were to tighten the tolerances you might need more input digits.
To clarify: I use solve() and vpasolve() here to show that there are solutions and to give us an idea of where they are so that we can explore what is needed in order to get fsolve() to work. solve() and vpasolve() are not intended to be part of the permanent solution (though if you have access to the Symbolic Toolbox, then solve() finds the exact solution easily)
It turns out that you need to be somewhat precise in order for fsolve() to be able to say it is happy, indicating that the equations are numerically sensitive.