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 }]

[![enter image description here][1]][1]

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.


  [1]: https://i.sstatic.net/nJ0OT.png