Supppose I have two states $$|\psi\rangle = \sum_{i=1}^n c_i |\psi_i\rangle\,,$$ $$|\phi\rangle = \sum_{i=1}^n c_i |\phi_i\rangle\,,$$ such that $|\psi_i\rangle$ and $|\phi_i\rangle$ are $\epsilon$ close in trace distance $\forall i$.
What is the upper bound on the trace distance between $|\psi\rangle$ and $|\phi\rangle$?
My intuition says it should be $n \epsilon$, but I can not prove it.
When taking your question literally, then no such bound can exist, i.e., there exist coefficients $\{c_i\}_i$ and element-wise trace-norm-close bases $\{\phi_i\}_i,\{\psi_i\}_i$ such that no constant $C>0$ satisfies $\|\,|\phi\rangle\langle\phi|-|\psi\rangle\langle\psi|\,\|_1\leq C\epsilon$. You even gave the counterexample yourself: $$c_1=c_2=\frac1{\sqrt 2},\quad \phi_1=\psi_1=|0\rangle,\quad \phi_2=|1\rangle,\quad \psi_2=-|1\rangle,\quad \epsilon=0\,.$$
[Original answer] If you relax the question and assume that $\psi_i,\phi_i$ are $\epsilon$-close not in trace norm, but rather in the usual Hilbert space norm, i.e. $\|\psi_i-\phi_i\|<\epsilon$, then one can prove that $$ \big\|\,|\psi\rangle\langle\psi|-|\phi\rangle\langle\phi|\,\big\|_1\leq 2\epsilon\sqrt n\tag1 $$ which is even better than the bound you conjectured (unless $n=2$ or $n=3$). To be clear, I do not claim that (1) is in any way optimal, but (1) is a nice bound because proving it is relatively simple. Said proof relies on the following observation: $$ \|\psi-\phi\|\leq\epsilon\sqrt n\tag 2 $$ This is most easily shown via the triangle inequality combined with Cauchy-Schwarz: \begin{align*} \|\psi-\phi\|&=\Big\|\sum_{i=1}^nc_i(\psi_i-\phi_i)\Big\|\\&\leq \sum_{i=1}^n|c_i|\,\|\psi_i-\phi_i\|\\&\leq\Big(\sum_{i=1}^n|c_i|^2\Big)^{1/2}\Big(\sum_{i=1}^n\|\psi_i-\phi_i\|^2\Big)^{1/2}\\ &\leq 1\cdot\Big(\sum_{i=1}^n\epsilon^2\Big)^{1/2}=\big(n\epsilon^2\big)^{1/2}=\epsilon\sqrt n\,. \end{align*} (here I assumed that $\sum_i|c_i|^2=1$ which is true, e.g., if the $\{\psi_i\}_i$ are orthonormal as is usually the case with these superposition expansions). An interesting side note is that this bound is optimal, that is, one cannot get rid of the $\sqrt n$ factor on the right-hand side.
In any case with this, (1) is just another application of the triangle inequality: \begin{align*} \big\|\,|\psi\rangle\langle\psi|-|\phi\rangle\langle\phi|\,\big\|_1&=\big\|\,|\psi\rangle\langle\psi|-|\psi\rangle\langle\phi|+|\psi\rangle\langle\phi|-|\phi\rangle\langle\phi|\,\big\|_1\\ &\leq\big\|\,|\psi\rangle\langle\psi|-|\psi\rangle\langle\phi|\,\big\|_1+\big\|\,|\psi\rangle\langle\phi|-|\phi\rangle\langle\phi|\,\big\|_1\\ &=\big\|\,|\psi\rangle\langle\psi-\phi|\,\big\|_1+\big\|\,|\psi-\phi\rangle\langle\phi|\,\big\|_1\\ &=2\|\psi-\phi\|\leq 2\epsilon\sqrt n\tag*{$\square$} \end{align*}
As a final remark, of course (1) only makes sense for $\epsilon$ small enough, i.e., it is useless for $\varepsilon\geq 1/\sqrt n$ as then, (1) becomes worse than the trivial bound $\|\,|\psi\rangle\langle\psi|-|\phi\rangle\langle\phi|\,\|_1\leq 2$.
