Actually, simply
p = First@NDSolve[{f''[x] == - Sum[z[i] cp[i], {i, 1, 2}],
f[10] == zf/(ep a^2), f'[20] == 0}, f, {x, 10, 100 }];
Plot[(f /. p)[x], {x, 10, 100 }]
gives the same answer. The term D[HeavisideTheta[x - 10], x] zf/(ep a^2) has no effect, because it is at the boundary. It is strange that NDSolve gives an error instead of just ignoring the term.