This problem can be solved by colocation method with using Euler wavelets. First, we map solution on unit interval $0\le x\le 1$. Then we define function f, first (f1) and second (f2) derivative as follows
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7; xmin = 0; xmax = 100; l = xmax - xmin; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[xmin + s*dx, {s, 0, nn}]; xcol =
Table[(xl[[s - 1]] + xl[[s]])/2, {s, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1];
intx1[y_] := Intx1 /. t1 -> y; intx2[y_] := Intx2 /. t1 -> y;
Psi[y_] := Psijk /. t1 -> y;
a[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
b[x_] := l^2 Piecewise[{{x^2, 0 <= x <= 1}, {0, True}}];
c1[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x1 = 1/l;
c2[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x2 = 10/l;
c3[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x3 = 5/l;
c4[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x4 = 15/l;
(*eq=D[a[x] f'[x]+b[x] f[x],x]+(c1[x+x1] f[x+x1]-c1[x] \
f[x])+(c2[x+x2] f[x+x2]-c2[x] f[x])+(c3[x-x3] f[x-x3]-c3[x] \
f[x])+(c4[x-x4] f[x-x4]-c4[x] f[x])\[Equal]0;bc1=f'[l]\[Equal]0;
bc2=f[l]\[Equal]0;*)
var = Array[v, {nn}]; f[x_] := var . intx2[x] + x v1 + v0;
f1[x_] := var . intx1[x] + v1; f2[x_] := var . Psi[x];
Finally we compute equation in collocation points and solve system of algebraic equations with FindRoot
eq = Table[
a[x] f2[x]/l^2 + f1[x]/l + b[x] f1[x]/l +
2 l x f[x] + (c1[x + x1] f[x + x1] -
c1[x] f[x]) + (c2[x + x2] f[x + x2] -
c2[x] f[x]) + (c3[x - x3] f[x - x3] -
c3[x] f[x]) + (c4[x - x4] f[x - x4] - c4[x] f[x]) == 0, {x,
xcol}];
eqs = Join[
Drop[eq, -1], {l^1.5 Sum[
Sqrt[xcol[[i]]] f[xcol[[i]]] (xl[[i + 1]] - xl[[i]]), {i,
nn}] == 1, f[1] == 0, f1[1] == 0}]; varM = Join[var, {v0, v1}];
sol = FindRoot[eqs, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]
Visualization f[x]
Plot[f[x/l] /. sol, {x, 0, l}, Frame -> True, PlotRange -> All]

Update 1. We also can use NMinimize
with bound f[x]>0
at 0<x<l
as follows (note we define exact integral intf[x]
instead of sum)
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 8; xmin = 0; xmax = 100; l = xmax - xmin; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[xmin + s*dx, {s, 0, nn}]; xcol =
Table[(xl[[s - 1]] + xl[[s]])/2, {s, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1]; Intx3 =
Integrate[Intx2 Sqrt[t1], t1, Assumptions -> t1 > 0];
intx1[y_] := Intx1 /. t1 -> y; intx2[y_] := Intx2 /. t1 -> y;
intx3[y_] := Intx3 /. t1 -> y; Psi[y_] := Psijk /. t1 -> y;
a[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
b[x_] := l^2 Piecewise[{{x^2, 0 <= x <= 1}, {0, True}}];
c1[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x1 = 1/l;
c2[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x2 = 10/l;
c3[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x3 = 5/l;
c4[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x4 = 15/l;
(*eq=D[a[x] f'[x]+b[x] f[x],x]+(c1[x+x1] f[x+x1]-c1[x] \
f[x])+(c2[x+x2] f[x+x2]-c2[x] f[x])+(c3[x-x3] f[x-x3]-c3[x] \
f[x])+(c4[x-x4] f[x-x4]-c4[x] f[x])\[Equal]0;bc1=f'[l]\[Equal]0;
bc2=f[l]\[Equal]0;*)
var = Array[v, {nn}];
intf[x_] := var . intx3[x] + 2/5 v1 x^(5/2) + 2/3 v0 x^(3/2);
f[x_] := var . intx2[x] + x v1 + v0; f1[x_] := var . intx1[x] + v1;
f2[x_] := var . Psi[x];
eq = Table[
a[x] f2[x]/l^2 + f1[x]/l + b[x] f1[x]/l +
2 l x f[x] + (c1[x + x1] f[x + x1] -
c1[x] f[x]) + (c2[x + x2] f[x + x2] -
c2[x] f[x]) + (c3[x - x3] f[x - x3] -
c3[x] f[x]) + (c4[x - x4] f[x - x4] - c4[x] f[x]), {x, xcol}];
bc =
Join[{l^(3/2) intf[1] == 1, f[1] == 0, f1[1] == 0},
Table[f[x] > 0, {x, xcol}]]; varM = Join[var, {v0, v1}];
solm =
NMinimize[{eq . eq, bc}, varM]
Visualization of solutions for nn=10, 16, 24, 28, 32, 56
We can combine all plots in one using Show

Update 2. We can use wavelets generated by $e^{-m t}$ to solve second problem as follows
UE[m_, t_] := Exp[-m t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[s*dx, {s, 0, nn}]; xcol =
Table[(xl[[s - 1]] + xl[[s]])/2, {s, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1 Sqrt[t1], t1, Assumptions -> t1 > 0];
intx1[y_] := Intx1 /. t1 -> y; intx2[y_] := Intx2 /. t1 -> y;
Psi[y_] := Psijk /. t1 -> y;
xmin = 0; xmax = 100; l = xmax - xmin; x0 = 1; x1 = x0/l;
dAfFun[x_] :=
Piecewise[{{1.4476805109594064` 10^-23 +
1.1735563659975857` 10^-24 x^(1.4), xmin <= x <= xmax}}, 0]
AfFun[x_] :=
Piecewise[{{4.17425813460343` 10^-24 + 4.4423429301477025` 10^-23 x,
xmin <= x <= xmax}}, 0]
Rx1Fun[x_] :=
Piecewise[{{4.537932520459608` 10^-17/(1 +
1.3188535755423771` Exp[-0.07679009872276484` (x - x0)]) -
4.529403364443408`*^-17/(1 +
1.3148412969171435` Exp[-0.07683234639209759` (x - x0)]),
x0 <= x <= xmax}}, 0]
a[x_] := dAfFun[x l] 10^20
b[x_] := AfFun[x l] 10^20
c1[x_] := Rx1Fun[x l] 10^20
var1 = Array[v, {nn}]; var2 = Array[u, {nn}];
intf[x_] := var1 . intx2[x] + 2/3 v0 x^(3/2);
f[x_] := var1 . intx1[x] + v0; f1[x_] := var1 . Psi[x];
g1[x_] := var2 . Psi[x]; g[x_] := var2 . intx1[x] + u0;
(*eq={a[x] f'[x]+b[x] f[x]==g[x],D[g[x],x]+(c1[x+x1] f[x+x1]-c1[x] \
f[x])==0};bc1=f'[l]\[Equal]0;
bc2=f[l]\[Equal]0;*)
eq =
Join[Table[g1[x]/l + (c1[x + x1] f[x + x1] - c1[x] f[x]), {x, xcol}],
Table[-g[x] + a[x] f1[x]/l + b[x] f[x], {x, xcol}]]; varM =
Join[var1, var2, {v0, u0}];
eqs = Join[
Table[eq[[i]] == 0, {i, 2, Length[eq] - 1}], {l^(3/2) intf[1] == 1,
f1[0] == 0, f[1] == 0, g[1] == 0}]; sol =
FindRoot[eqs, Table[{varM[[i]], 1/10}, {i, Length[varM]}]];
Visualization numerical solution and PDF[ChiDistribution[1], x]
{Plot[PDF[ChiDistribution[1], x], {x, 0, 10}, PlotRange -> All,
Frame -> True, PlotLabel -> "PDF[ChiDistribution[1],x]"],
Plot[Evaluate[{f[x/l] /. sol}], {x, 0, 20}, Frame -> True,
PlotRange -> All, PlotLabel -> nn],
Plot[Evaluate[{f[x/l] /. sol}], {x, 0, 100}, Frame -> True,
PlotRange -> All, PlotLabel -> nn]}

Note, that we can solve linear system of equations with LinearSolve
using
{vec, mat} = CoefficientArrays[eqs, varM];
soll = LinearSolve[mat, -vec];
rule = Table[varM[[i]] -> soll[[i]], {i, Length[varM]}];
Visualization two numerical solutions
Plot[Evaluate[{f[x/l] /. rule, f[x/l] /. sol}], {x, 0, l},
Frame -> True, PlotRange -> All, PlotLabel -> nn,
PlotLegends -> {"LinearSolve", "FindRoot"},
PlotStyle -> {Blue, Dashed}]

bc1 = f'[l] == 0; bc2 = f[l] == 0
solution is $f(x)=0$. $\endgroup$