Can the Helmholtz equation be applied for a system with multiple frequencies?

If I understand what you want to do correctly, then yes, you can do it; specifically because the wave equation is linear. If you have two solutions to the wave equation $$ psi_1(vec{r},t) = phi_1(vec{r}) e^{i omega_1 t}, qquad psi_2(vec{r},t) = phi_2(vec{r}) e^{i omega_2 t}, $$ then the sum of these two solutions $psi = psi_1 + psi_2$ is also a solution to the wave equation: $$ nabla^2 psi + frac{1}{v^2} frac{partial^2 psi}{partial t^2} = 0. $$ But $phi_1$ and $phi_2$ will each satisfy the Helmholtz equation for a particular value of $k$: $$ nabla^2 phi_1 = k_1^2 phi_1, qquad nabla^2 phi_2 = k_2^2 phi_2 $$ with $k_i = omega_i/v$.

The full solution will then be $$ psi(vec{r}, t) = phi_1(vec{r}) e^{i omega_1 t} + phi_2(vec{r}) e^{i omega_2 t}. $$ Note, however, that these two waves don't "interfere" with each other per se, since they're at different frequencies. For example, the time-averaged amplitude of the wave at a point $vec{r}$ will be $frac{1}{2} (|phi_1|^2 + |phi_2|^2)$, rather than $frac{1}{2} |phi_1 + phi_2|^2$.

In brief: if you're looking at solutions more complex that a standing wave, you may have to resort to the time-dependent wave equation. However, we can use to our advantage that any general solution can be expressed as the weighted sum of all mode shapes, where the weights are functions of time.

I will assume that your problem is an initial conditions problem, where you know what the displacement and velocity of the system is at $t=0$, and you want to find how the system evolves over time. If not, please let me know.

Solving for propagating waves

Note that the Helmholtz equation

$$nabla^2 u = -lambda u$$

will have multiple (if not infinite) solutions, where each solution corresponds to each standing wave. Each solution consists of a pair made up of an eigenfunction/mode shape $u_n(mathbf{x})$ and an eigenvalue $lambda_n$, where $n$ denotes the $n^mathrm{th}$ standing wave. The eigenvalue $lambda_n$ is related to resonant frequency $omega_n$ by

$$lambda_n = frac{omega_n^2}{c^2}$$

If we want to understand how waves other than standing waves - such as propagating waves - behave, we do need to resort to using the time-dependent wave equation again:

$$ddot u = c^2 nabla^2 u$$

Trying to numerical solve the wave equation from scratch is very cumbersome, involving solving over both space and time domains.

However, we can greatly simplify things if we know all the mode shapes and resonant frequencies. A general solution of the wave equation can be expressed as the weighted sum of the mode shapes:

$$u(mathbf{x},t) = sum_{n=1}^{N} q_n(t) , u_n(mathbf{x})$$

$N$ is the number of modes, which will be infinite for continuous systems like strings, water waves, etc.

Now, to solve for the general solution, we just need to determine these "weighting functions"/modal coordinates $q_n(t)$. Note that the case of a standing wave is when all but one of the modal coordinates are zero.

Provided that the mode shapes are normalised, say, such that

$$int_Omega |u_n(mathbf{x})|^2 , mathrm{d}m = 1$$

then we can transform our more complicated time-and-space-dependent wave equation into individual simpler time-dependent ODEs:

$$ddot{q}_n + omega_n^2 q_n = 0$$

This is the ODE for undamped harmonic motion, at frequency $omega_n$.

There will exist an ODE for each modal coordinate, and these are much simpler to numerical solve than a PDE over space and time. In fact, for this simple case, analytical expression can be obtained for the general solution of each modal coordinate, and you would only have to solve for the constants of the general solution.

We need the values $q_n(0)$ and $dot{q}_n(0)$ (the initial conditions) of each ODE to be able to solve, and these can be determined from $u(x,0)$ and $dot{u}(x,0)$ using the following expressions by setting $t=0$:

$$q_n(t) = int_Omega u(mathbf{x},t) , u_n(mathbf{x}) , mathrm{d}m$$

$$dot{q}_n(t) = int_Omega dot{u}(mathbf{x},t) , u_n(mathbf{x}) , mathrm{d}m$$

Depending on whether the mode shapes are known as analytical or numerical expression, these expressions may need to be evaluated by analytical or numerical integration.

In summary

In summary, to find the solution for propagating waves:

• Find ALL the mode shapes and resonant frequencies using Helmholtz equation
• Determine the boundary conditions for all the ODEs for each mode
• Solve each ODE for $q_n(t)$
• Construction the general solution using $u(mathbf{x},t) = sum_{n=1}^{N} q_n(t) , u_n(mathbf{x})$

Big caveat

Note that we require knowledge of ALL mode shapes and resonant frequencies. For continuous system like water waves, there are infinitely many modes, and so this task is infeasible. However, a practical solution is to approximate using finitely many modes by excluding modes whose resonant frequency are well above the frequency range of interest, effectively applying a low pass filter to the true solution.