What you are asking about falls into the realm of spectral estimation, of which frequency estimation is a particular case.
In general, with no prior knowledge about the frequencies in your original signals, you cannot do what you ask. One fundamental reason is aliasing. If you sample at some frequency \$f_s\$, you cannot distinguish input frequencies \$nf_s - f_0\$, \$nf_s + f_0\$ for different n from 0 to infinity. For example, if you sample at 100 Hz, you can't tell the difference between a 20 Hz input, an 80 Hz input, 120 Hz, 180, 220, etc.
Another limitation, if you know nothing about your input frequencies, is that the precision with which you can estimate your input frequencies is limited by the length of time you sample for. For example, if you sample at say 100 Hz for 1 s, you might (very roughly speaking) just barely be able to distinguish between a 20 Hz input and a 21 Hz input. If you sample for 100 s, you might barely be able to distinguish between a 20 Hz input and a 20.01 Hz input, etc.
These limitations apply whether you have an input that is formed by combining two independent sources or if you just have a single, pure, sinewave input and you want to estimate its frequency. Of course the second limit has some relevance if you have two inputs at closely spaced frequencies and you want to be able to separate them.
You may also be able to get some help with this at dsp.stackexchange.com, though the typical answer over there requires a substantial amount of mathematical background that may make them difficult to understand.
This is a windowing artifact.
The linked code pads out a 10,000 sample signal with zeroes so that the length is a power of two.
%% Author :- Embedded Laboratory
%%This Project shows how to apply FFT on a signal and its physical
% significance.
fSampling = 10000; %Sampling Frequency
tSampling = 1/fSampling; %Sampling Time
L = 10000; %Length of Signal
t = (0:L-1)*tSampling; %Time Vector
F = 100; %Frequency of Signal
%% Signal Without Noise
xsig = sin(2*pi*F*t);
...
%%Frequency Transform of above Signal
subplot(2,1,2)
NFFT = 2^nextpow2(L);
Xsig = fft(xsig,NFFT)/L;
...
Note that in the above code, the FFT is taken with the FFT size NFFT which is the next power of 2 bigger than the signal length (in this case, 16,384.) From the Mathworks fft() documentation:
Y = fft(X,n) returns the n-point DFT. fft(X) is equivalent to fft(X, n) where n is the size of X in the first nonsingleton dimension. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in the same manner.
This means that you are not actually taking a FFT of a 'pure sine wave' - you are taking the FFT of a sine wave with a flat signal after it.
This is equivalent to taking the FFT of a sine wave multiplied with a square window function. The FFT spectrum is then the convolution of the sine wave frequency spectrum (an impulse function) with the square wave frequency spectrum (sinc(f).)
If you change L = 16,384 so that there is no zero-padding of the signal, you will observe a perfect FFT.
Further search keywords: "Spectral Leakage", "Window Function", "Hamming Window".
Edit: I cleaned up some material I wrote on this topic back in university, which goes into substantially more detail. I've posted that on my blog.
Best Answer
The Fourier transform is linear, which means that if you sum two signals, then also their spectrum will be the sum.
$$ \mathscr{F} [x(t) + y(t)] = \mathscr{F}[x(t)] + \mathscr{F}[y(t)] $$
If you consider the power, you can have two cases:
The signals have different frequency
Then their power just sums up, because: $$ P_{tot} = P_1 cos (\omega _1 t + \phi _1) + P_2 cos (\omega _2 t + \phi _2) $$
The signals have the same frequency
Then if they are in phase you have: $$ P_{tot} = P_1 cos (\omega _0 t + \phi _0) + P_2 cos (\omega _0 t + \phi _0) = (P_1 + P_2) \cdot cos(\omega _0 t + \phi_0) $$ If they're out of phase, the resulting power will be lower and precisely: $$ P_{tot} = P_1 cos (\omega _0 t + \phi _1) + P_2 cos (\omega _0 t + \phi _2) $$ and setting arbitrarily \$\phi_1 = 0\$ we obtain: $$ P_{tot} = P_1 cos (\omega _0 t) + P_2 cos (\omega _0 t + \Delta \phi) $$ which for \$\Delta \phi = \pi\$ gives the subtraction of the signals.