Consider Schroedinger equation, which I write in the form $$ (\mathscr{L}+V)\psi=0$$where $\mathscr{L}$ is the kinetic and time-derivative operator. Now, imagine I have two point sources 1,2 with potential $V_i(r)=f_i(r) \delta (r-r_i)$. Let $\psi_1$ (resp $\psi_2$) be the solution when I switch on only source 1 (resp 2). It is clear that when both sources are switched on then potential $V=V_1+V_2$. However, unlike the sum of amplitudes, $\psi\neq \psi_1+\psi_2$ , instead it is $$\psi=\psi_1+\psi_2-\frac{1}{\mathscr{L}+V}(V_1\psi_2+V_2\psi_1)$$ Consequently, every other quantities that can be derived from $\psi$ won't respect the usual linearity. I could use similar setup in double-slit experiment , where $V_i$ is infinite except at the location of slit $i$ (The total potential $V$ is infinite except at location of slits 1 and 2). The argument above equally applies for classical theory like Electro-magnetism (as shown here).
It seems from above argument that if EM field (or wavefunction) from multiple sources are given, I cannot always vectorially add them to get the overall field strength at a point. So when should the vector sum of field strength/ sum of amplitudes be valid?