This is more like a comment.
You have a discontinuous ode, so to turn off the discontinuity during NDsolve's processing of your ode with DiscontinuityProcessing,
NDSolve[{f''[x]==DiracDelta[-10 + x] zf/(ep a^2) - 1 Sum[z[i] cp[i], {i, 1, 2}],
f[10] == zf/(ep a^2), f'[20] == 0}, f, {x, 10, 100},
Method -> {"DiscontinuityProcessing" -> False}]
which generate this error,
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 10.`.
So I added a submethod Method -> "ExplicitEuler"
p = NDSolve[{f''[x] ==
DiracDelta[-10 + x] zf/(ep a^2) - 1 Sum[z[i] cp[i], {i, 1, 2}],
f[10] == zf/(ep a^2), f'[20] == 0}, f, {x, 10, 100},
Method -> {"FixedStep", Method -> "ExplicitEuler",
"DiscontinuityProcessing" -> False}]
which produced a solution but with a warning,
NDSolve::nlnum: The function value {-0.0000234567,2.5*10^-6+0.00125 DiracDelta[0.]} is not a list of numbers with dimensions {2} at {x,f[x],(f^[Prime])[x]} = {10.,0.00125,-0.0000234567}.
Plot[f[x] /. p, {x, 10, 100}]
