You made me curious with this example.
(Actually it's obvious that you started from an approximation of ln(x)  [maybe using remez].)

So, you want the roots of f  in the interval (0,1]

It is possible to show that f has exactly 35 roots in 0..1 and they can be computed,
but unfortunately the Maple built-in commands are useless for this example. 
Especially RootFinding:-NextZero  gave very poor results.
Disappointing!

restart;
f := x -> 42312393/170170 + (87300630621*x)/5005 + (84260074354272*x^2)/1001 + (11572751542512000*x^3)/143 + (311217957498451500*x^4)/13 + (13736312974717541208*x^5)/5 + (702361109129611835904*x^6)/5 + (24407857262955082295808*x^7)/7 + (309128866743376578380625*x^8)/7 + (2034476230680567168673750*x^9)/7 + 971145727536841731347616*x^10 + (15981733578778631623709568*x^11)/11 + (3759286855176800298957060*x^12)/11 - (191464620989481690770115000*x^13)/143 - (1305974159036375354989560000*x^14)/1001 - (37785618862730180432031744*x^15)/91 - (8073200612643819138676179*x^16)/182 - (231388810612770205973655*x^17)/221 + (190600129650794094000*x^17 + 13377406242490734054600*x^16 + 195343515041861129011200*x^15 + 1040537189498550048000000*x^14 + 2518477612889131719000000*x^13 + 3093773724805440474252000*x^12 + 2048291569526360589849600*x^11 + 753720641133895782547200*x^10 + 155778902350755576750000*x^9 + 17974488732779489625000*x^8 + 1132669320760996761600*x^7 + 37459215415250044416*x^6 + 610748077422555072*x^5 + 4463801814780000*x^4 + 12619635840000*x^3 + 10816830720*x^2 + 1779084*x + 18)*ln(x):
plot(f, 1e-8.. 2*10^(-6))

isolate(f(x), ln(x));
g:=unapply((lhs-rhs)(%), x);

g1:=normal(diff(g(x),x)):
sturm(numer(g1),x,0,1),  sturm(denom(g1),x,0,1);

                             34, 0

# So, g' has 34 real roots and all are in (0,1);
 

Digits:=100;
L:= [10.^(-6), fsolve(numer(g1)),  1.];
signum~(g~(L));nops(%);
# [-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1,  -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, #  , -1, 1]
#                               36
# So, g (and f) has exectly 35 roots in (0,1);

for i to 35 do i=fsolve(g, L[i]..L[i+1]) od;

1 = .191561287749093307937227033160539875133752143080898731095086117213834851920515705439e-5
2 = .386863114821857362406592315527508687579950813861777293116197038157920498833807772632e-4
3 = .224368589492033177901033815331319924131620771048548665891126694787181114809420099728e-3
4 = .757628876465743496019898819891689451187953460530807054577709360768356026307038622925e-3
5 = .191586043560551126248587315094560734538794188822600052918532582221573212092426725234e-2
6 = .404341050752997444506457044613767563622770741756842016619777344000424267695214131236e-2
7 = .753803071564053290619849328637704614199820872264894616553942256990815977878634237764e-2
8 = .128347814103994790541831995624033323707358509644963097185620852545664106970618325968e-1
9 = .203881744484780444436159799720467561935237060222424328004944514734354971452828369636e-1
10 = .306531122499106815328715860139756718990911721033357327791034227500491192932066198881e-1
11 = .440652608872939466182326564099001428304047894621422939893923731352872694517742934307e-1
12 = .610215129442692726469217238624703315785764842073204122632174231202851842298289275127e-1
13 = .818611973531043298152537221149426167590700088952157881607430662062960493383888023674e-1
14 = 106848673577559840945862773565531573854056603799252847848710761297768054444573771543
15 = 136157907371746382188673002171784836243338456425810338155943739653948751569648355407
16 = 169859566116596553194242901269900419712969903139816354118873867296351077315658311599
17 = 207911095254057736286977980962266715752752350283085803618405408714970551082808317735
18 = 250150145986443098687961663490645213960126759514191210023876415880465673699179966025
19 = 296291621034620997624748934217305561232655042486558270046438747007382649214020765990
20 = 345928493082594386040296659177343956611809301967460342049504180924850028932218436257
21 = 398536433083337831429021148180533317041041457085873035259687560041254081758343858771
22 = 453482166534225841822476226178659839805335530836077276359472418786766336138346293288
23 = 510035358874306638935680974371652097595533372440279658507373385287879407012832470726
24 = 567383719969290157155679436652208052931514611589201379988400501404970528045712346664
25 = 624650915713400311490990374296520507166356685922856856979075726522784716531016773873
26 = 680916785283116652996592071240106370085097908173787711070599706365125792126746060999
27 = 735239288335241089539005060331848708632199246169167693175129656187672273735745985465
28 = 786677549795751681675503376421685045544117797715570060193731242723710977234911696913
29 = 834315332668459220663475591628563704363467776558796639367542621477772723238905506446
30 = 877284252847918207927215084652553361226192179200405713639782815744067137857243235030
31 = 914786055398803215378013441256970465033525746848992402044164144886322908035392340196
32 = 946113301720552053182303564465904800429669666884353849407725334749814613284598367425
33 = 970667885708695016778657128606495810300922043980656327177238682949045151388078173429
34 = 987977026592788760481681524535038678110293678929974123268967108570856367223442229329
35 = 997708817965433194699497427026838774265049425782232856894868033444326438291042212089

r[0]:=1e-6;
for i do
r[i]:=RootFinding:-NextZero(g,r[i-1]+1e-10, maxdistance=1);
until r[i]=FAIL;      # disappointing!

# aborted !

Edit.  However, using 
RootFinding:-NextZero(g,r[i-1]+1e-10, maxdistance=1, guardDigits=30);
NextZero works!

Note that the parametrization is 4*Pi - periodic in u, so, no need to plot over 0..6*Pi.
 

restart;
f:=[arctan(-tan(u/4)) + ((cosh(v) - sinh(v))*sin(u))/2, -v/2 + ((cosh(v) - sinh(v))*cos(u))/2, 2*(cosh(v/2) - sinh(v/2))*sin(u/2)]:
e:=1.e-3:
P1:=plot3d(f, u = 0 .. 2*Pi-e,    v = -1 .. 1, grid = [80, 15], orientation = [50, 79], shading = XYZ, style = patch):
P2:=plot3d(f, u = 2*Pi+e .. 4*Pi, v = -1 .. 1, grid = [80, 15], orientation = [50, 79], shading = XYZ, style = patch):
plots:-display(P1,P2);

I simply insert (when needed) some empty execution groups.

discont(Re(p1), alpha): evalf(%)[]:
A:=0, %[-3..-1], 5:
plots:-display(seq(plot([Re(p1),Re(p2)],alpha=A[i]+0.001 .. A[i+1]-0.001, view=-2..2, color=[red,blue]), i=1..4));

Download dsolveD.mw

It seems that you made yourself (using a palette?) a variable named alpha~. Don't do that, just use alpha. Actually it is much better to use % (or the label of the expression) rather than copy+paste or re-type the expression.
simplify(fracdiff(%, t, alpha));

Compute the series around x = 0. You will see that Maple's result is correct.

The procedure is very simple and practically it cannot be made faster (the bottleneck being isprime).
For n>28 you will need a lot of patience. But you can save f and rerun it another day, the computed values are saved too.
 

restart;
f:=proc(n::posint)
local a,b,c;
a,b:=f(n-2),f(n-1);
while not isprime((c:=a+b)) do a,b:=b,c od;
if n<13 then lprint(n=c) else
   printf("%d = %s ... %s len=%d\n", n, ""||c[1..10], ""||c[-10..-1], length(c)) fi; 
f(n):=c;
end proc:
f(1):=2:f(2):=3:

f(27):

3 = 5
4 = 13
5 = 31
6 = 313
7 = 2659
8 = 96979
9 = 97340263
10 = 96133996771
11 = 288596670839
12 = 35613385860024917251
13 = 1210855125 ... 1377274153 len=22
14 = 4191695536 ... 7350583473 len=23
15 = 1559140836 ... 5674347327 len=35
16 = 1169745412 ... 9762684239 len=40
17 = 3996415043 ... 1378463323 len=55
18 = 4535544101 ... 7757493097 len=64
19 = 2625668239 ... 2056367381 len=113
20 = 2276456306 ... 0210477381 len=138
21 = 2046360541 ... 3641285093 len=168
22 = 4162619505 ... 8047946001 len=271
23 = 1930123412 ... 5467084469 len=276
24 = 1665092783 ... 3398723741 len=287
25 = 7550383970 ... 9920377053 len=421
26 = 1325989112 ... 9358807409 len=691
27 = 1475647825 ... 1633573333 len=1030
 

A = A(i) can be proved to be equal to - (-d)^i. Use induction. Maple does not want to simplify A(i) to -(-d)^i. It would be possible to use Maple for A(i+1) etc, but it's easier by hand..

The rest is simple, because B = - A(i+1)/d = - (-d)^i = A(i).

A:= t -> numapprox:-infnorm((f_exact-f_numeric)(x,t), x=0..2):
numapprox:-infnorm('A'(t), t=0..2);
#                         0.02681084520

AA:=convert(A,rational):
bb:=convert(b,rational);
xx := LinearSolve(AA, bb);

restart;
g[2]:=0: g[1]:=1: g[3]:=1: alpha:=2: c:=1:
k:=(g[2]+2*g[3]-3*g[1]*alpha)/(6*g[1]*g[3]):
omega:=(((1-3*g[1]*k)*(2*k-c-3*g[1]*(k^2))  )/(g[1]))+(k^2)-g[1]*(k^3):
uu11:=1/(g[2]+2*g[3])^(1/2)*(-3*(3*k^2*g[1]+c-2*k))^(1/2
)*sin(1/2/g[1]*2^(1/2)*(g[1]*(3*k^2*g[1]+c-2*k))^(1/2)*(-c*t+x))/
cos(1/2/g[1]*2^(1/2)*(g[1]*(3*k^2*g[1]+c-2*k))^(1/2)*(-c*t+x))*exp(
I*(k*x-omega*t)):
pde:=I*Diff(u(x,t),t)+Diff(u(x,t),x$2)+alpha*(abs(u(x,t))^2)*u(x,t)+ 
I*( g[1]*Diff(u(x,t),x$3) + g[2]*(abs(u(x,t))^2)*u(x,t) + g[3]*Diff((abs(u(x,t))^2),x)*u(x,t) ):
eval(pde, u(x,t)=uu11):
Z:=value(%):
eval(Z, [x=1,t=2]): evalf(%);  # <>0
eval(Z, [x=2,t=1]): evalf(%);  # <>0

                  -12687.93889 - 28829.95560 I

                  25022.56665 - 19131.68293 I

The two functions are equal (f and the spline). This is normal, f being a polynomial of degree <=3). So, if you want to test your construction, choose another f.
convert(NaturalSpline(x) - f, rational) ; # ==> 0.

BTW. Your f is an expression, not a procedure. So, don't use f(x) because it's a nonsense. 

So, you probably want to approximate a solution of ODE0 by U0. Why not something like this?

restart;
ODE0:= diff(U(x),x$2) + A*U(x) + B*U(x)^3;
U0:=x -> sin(mu*x)/(K + L*cos(mu*x));
Z:=eval(ODE0, U=U0);
series(Z,x,4);
coeffs(convert(%,polynom),x);
solve([%],[K,L],explicit);