I got general solution for PDE I solved
$$f(x, y) = \sum_{n=0}^{\infty} \exp \left(-A \frac{n^2 \pi^2}{4L}y\right)\left[ B_n \cos\left(\frac{2n\pi}{2L}x \right) + C_n\sin\left(\frac{(2n+1)\pi}{2L}x \right)\right] $$
with boundary contition
$$f(x, 0) = N\delta(x-x_0)$$
where $A \in \mathbb{R}_+$, $n \in \mathbb{Z}$, and $B_n$ and $C_n$ are constants related to $n$, $2L$ is length of '$x$' interval where $f(x, y)$ is defined. $\delta(x, x_0)$ is Dirac's Delta
$x \in [-L; L]$
$y \in [0; \infty)$
I'm hoping to designate $B_n$ and $C_n$, but as far as I know I can not get values of two constants from one boundary condition.
What should I do? Can I limit my sum only to odd/even numbers (then either 'sin' term, or 'cos' will survive) without losing generality of my solution?
If so, how should I proceed then?
POSSIBLE SOLUTION Setting $x_0 = 0$ I integrate both sides of $$f(x, 0) = \sum_{n=0}^{\infty} \left[ B_n \cos\left(\frac{2n\pi}{2L}x \right) + C_n\sin\left(\frac{(2n+1)\pi}{2L}x \right)\right]$$ with $$\int dx \cos\left(\frac{2n\pi}{2L}x \right)$$ Using Delta proprety LHS is $$N \int dx \delta(x) \cos\left(\frac{2n\pi}{2L}x \right) = N$$ RHS using orthogonality of trigonometric functions(since $2n+1$ isn't equal to $2m$) will be $$\sum_{n=0}^{\infty} \frac{B_n L}{2 \pi}\delta_{nm} $$ where $\delta_nm$ is Kronecker's delta. Summing with delta leads to $$N = \frac{B_m L}{2 \pi} \rightarrow B_m = N \frac{2\pi}{L}$$
Analogously integrating with $$\int dx \sin\left(\frac{(2m+1)\pi}{2L}x \right)$$ will lead to $$0 = \frac{C_m L}{2 \pi} \rightarrow C_m = 0$$
Then, final soultion will be $$f(x, y) = \sum_{n=0}^{\infty} \exp \left(-A \frac{n^2 \pi^2}{4L}y\right)\left[ N \frac{2\pi}{L} \cos\left(\frac{2n\pi}{2L}x \right)\right] $$
I assume you meant to write $f$ instead of $f_n$. The $n$ is a summation index.
At $y=0$ your function is a Fourier series written in terms of cosines and sines. What you need is the Fourier series of a delta function, and you can then match the coefficients. The Fourier series of the standard delta function centered at the origin is easily found online or in books, and you can then use the standard formula to shift the function. Pay attention to different normalizations of the Fourier series, though.
Alternatively, you can compute that $$ \int_{-L}^Lf(x,0)\cos(2m\pi x/2L)dx = \frac{LB_m}{2\pi} $$ and similarly for sines. If you use this formula for $f(x,0)=N\delta(x-x_0)$, you get formulas for $B_n$ and $C_n$.
To get the alternate formula, compute the integral term by term (commuting the integral and the sum) and compute the integrals of $\cos(2m\pi x/2L)\cos(2n\pi x/2L)$ (and similarly for sines) by hand. You will see that almost all the integrals are zero. Your cosines and sines have slightly different arguments, but that is not an issue. This is at the core of Fourier series, and I suggest reading on the topic if you are not familiar. Figuring out the integrals does take a while, and this is part of the work of building up a theory of Fourier series.