Free particle and infinite well, why do solutions proceed differently? When one solves the Schrodinger's equation (I completely neglect the time part of wave function in all text) for a free particle, one ends up with fundamental system (for space variable):
$$FS=\{e^{ikx},e^{-ikx}\}$$ and thus one should write the general solution for the spatial part of wave function as:
$$\psi(x)=Ae^{ikx}+Be^{-ikx}$$
I have not seen this done anywhere. Every text ends up with fundamental system, and takes the individual exponentials as particular solutions and recasts it:
$$\psi(x)=e^{\pm i k x} = e^{ipx/\hbar}\,,$$
where $p\equiv\pm \hbar k$. I should mention even before starting out, the texts like to say as energy is non-negative, thus they can denote $k^2=\frac{\hbar^2}{2m}E$ and then they claim that $k$ is a real non-negative number (why? I have no idea)
To construct the general solution they say, that any solution can be a linear combination of the $e^{ipx/\hbar}$. Ie:
$$\psi(x)=\frac{1}{\sqrt{2\pi\hbar}}\int_{-\infty}^{+\infty} c(p) e^{ipx/\hbar}\,{\rm d}p$$
(yes the integrand should have the time part here, but this is irrelevant for my mental problem).
On the other hand, when the text go on to solve infinite potential well, they at first glance do something completely different, despite the problem, within the well, is the very same, ie. within the well, one solves the very same differential equation which also gives the fundamental system as above. But in this case, they say, that the general solution is :
$$\psi(x)=Ae^{ikx}+Be^{-ikx}$$
and then they go on to show quantization arising from boundary conditions (but this is no longer interesting to me now).
What I cannot seem to be able to understand is the seemingly different approach to the very same problem (again quantization is something that happens after one says that one takes such and such general solution).
I am not stating that something is wrong, I just cannot reconcile the different approaches, ie. why in the case of free particle does one takes only the particular solutions (the fact that they are eigenfunctions of impuls operator and the general solution is not, is not a good reason, I think, just a convenient consequence) and in the case of infinite well, the solution proceeds as it (from my perspective) should, getting a general solution.
 A: This question is fundamentally about how to label the solutions to eigenvalue problems, and it's largely a notation issue.$^{[a]}$
Eigenvalue problem
The generic eigenvalue problem is
$$T \left \lvert \lambda \right \rangle = \lambda \left \lvert \lambda \right \rangle$$
where $T$ is a linear transformation, $\lambda$ is the eigenvalue, and $\left \lvert \lambda \right \rangle$ is the eigenvector.
When working with a particular linear transformation, it's very common to label vectors by their eigenvalues under that transformation.
Suppose we find two solutions $\left \lvert \lambda_1 \right \rangle$ and $\left \lvert \lambda_2 \right \rangle$ that both satisfy the equation and with the same eigenvalue, i.e.
$$T \left \lvert \lambda_1 \right \rangle = \lambda \left \lvert \lambda_1 \right \rangle \quad \text{and} \quad T \left \lvert \lambda_2 \right \rangle = \lambda \left \lvert \lambda_2 \right \rangle \, .$$
Because $T$ is linear, any weighted sum of $\left \lvert \lambda_1 \right \rangle$ and $\left \lvert \lambda_2 \right \rangle$ is also a solution:
$$
T\left( a \left \lvert \lambda_1 \right \rangle + b \left \lvert \lambda_2 \right \rangle\right)
= a T \left \lvert \lambda_1 \right \rangle + b T \left \lvert \lambda_2 \right \rangle
= a \lambda \left \lvert \lambda_1 \right \rangle + b \lambda \left \lvert \lambda_2 \right \rangle
= \lambda (a \left \lvert \lambda_1 \right \rangle + b \left \lvert \lambda_2 \right \rangle)\, . $$
Free particle
In this case, we want to solve the eigenvalue problem
$$\left( -\frac{\hbar^2}{2m} D_x^2\right) \left \lvert \psi \right \rangle = E \left \lvert \psi \right \rangle \tag{$\star$}$$
where $D_x$ means "derivative with respect to $x$.
We can rewrite this equation in the $x$ basis if we want, by taking the inner product of both sides with $\left \langle x \right \rvert$:
\begin{align}
\left \langle x \right \rvert \left( -\frac{\hbar^2}{2m} D_x^2\right) \left \lvert \psi \right \rangle
&= \left \langle x \right \rvert E \left \lvert \psi \right \rangle \\
\left( - \frac{\hbar^2}{2m} \right)\frac{d^2}{dx^2} \psi(x) &= E \psi(x)
\end{align}
where we defined $\psi(x) \equiv \left \langle x \right \rvert \psi \rangle$.
This is the usual form of the Schrodinger equation, written in the position basis.
It's clear that this is a linear transformation because the derivative itself is linear, i.e. $D_x^2(f + g) = D_x^2 f + D_x^2 g$.
There are two solutions to this equation:
$$\exp \left(i \frac{\sqrt{2mE}}{\hbar} x \right)
\quad \text{and} \quad
\exp \left(-i \frac{\sqrt{2mE}}{\hbar} x \right) \, .$$
It's convenient to define $k \equiv \sqrt{2mE}/\hbar$ and write the two solutions as $\exp \left( \pm i k x \right)$.
Note that these vector solutions to equation $(\star)$ are written in the $x$ basis.
We could more generally write them as $\left \lvert \pm k \right \rangle$, in which case $ \left \langle x | k \right \rangle = \exp(i k x) \equiv \psi_k(x) \, .$

Spectral theorem
It is also the case that the eigenvalues of a Hermitian linear transformation form a basis for the vector space.
This is called the spectral theorem.

The linear operator $D_x^2$ is Hermitian (either take my word or check it yourself!).
Therefore, from the spectral theorem, we know that can write any quantum state as a weighted sum of $\left \lvert k \right \rangle$ states:
\begin{align}
\left \lvert \psi \right \rangle
&= \int_{-\infty}^\infty \frac{dk}{2\pi} \underbrace{\tilde \psi(k)}_\text{weights} \left \lvert k \right \rangle \\
\langle x \left \lvert \psi \right \rangle
&= \int_{-\infty}^\infty \frac{dk}{2\pi} \tilde \psi(k) \left \langle x \lvert k \right \rangle \\
\psi(x)
&= \int_{-\infty}^\infty \frac{dk}{2\pi} \tilde \psi(k) e^{ikx} \, . \tag{$\star \star$}
\end{align}
The sum goes over negative and positive values of $k$ because both $\left \lvert k \right \rangle$ and $\left \lvert -k \right \rangle$ are distinct, linearly independent eigenvectors in the basis set.
If we want a wave function $\left \lvert \psi \right \rangle_E$ that has a particular value of $E$, then we'd restrict to just those values of $k$ that match the desired value of $E$, i.e. defining a new symbol $K = \sqrt{2mE}/\hbar > 0$, we'd have
$$\left \lvert \psi \right \rangle_E = \tilde \psi(K) \left \lvert K \right \rangle + \tilde \psi(-K) \left \lvert -K \right \rangle \, . \tag{$\star \star \star$}$$
Now let's qnswer your questions

and thus one should write the general solution for the spatial part of wave function as:
  $$ \psi(x)=A e^{ikx} + Be^{−ikx}$$

Yes, that's the general solutions.
If $k = \sqrt{2mE}/\hbar$, then what you've written there is the general solution to the free particle Schrodinger equation with eigenvalue $E$.

I have not seen this done anywhere.

Well, at least you've seen it in Equation $(\star \star \star)$.
Start with $(\star \star \star)$, define $A\equiv \tilde \psi(K)$, $B \equiv \tilde \psi(-K)$, change notation $K \to k$, and take the inner product of both sides with $\left \langle x \right \rvert$ to get exactly the same equation as you expected.

Every text ends up with fundamental system, and takes the individual exponentials as particular solutions and recasts it:
  $$ \psi (x)= e^{\pm ikx} = e^{ipx/\hbar} \, ,$$
  where $p\equiv \pm \hbar k$.

That's the textbook author being sloppy and imprecise.
In this case, the author is using the symbol $\psi(x)$ to indicate "one single term in the expansion of the wave function $(\star \star)$ for a given energy $E=\hbar^2 k^2 / 2m$".
The author may be doing this because they are thinking of a particle moving in a particular direction for some physical reason, e.g. perhaps they are thinking of a particle wave being emitted from a source.

I should mention even before starting out, the texts like to say as energy is non-negative, thus they can denote $k^2=\hbar^2 2mE$ and then they claim that $k$ is a real non-negative number (why? I have no idea)

This is, again, a notational sloppiness.
In $(\star \star)$, we have a sum over all positive and negative $k$.
However, we can rewrite $(\star \star)$ like this:
\begin{align}
\psi(x)
&= \int_{-\infty}^\infty \frac{dk}{2\pi} \tilde \psi(k) e^{ikx} \\
&= \int_{-\infty}^0 \frac{dk}{2\pi} \tilde \psi(k) e^{ikx}
 + \int_0^\infty \frac{dk}{2\pi} \tilde \psi(k) e^{ikx} \\
&= \int_0^\infty \frac{dk}{2\pi} \tilde \psi(-k) e^{-ikx}
 + \int_0^\infty \frac{dk}{2\pi} \tilde \psi(k) e^{ikx} \\
&= \int_0^\infty \frac{dk}{2\pi} \left( \tilde \psi(k) e^{ikx} + \tilde \psi(-k) e^{-ikx} \right) \, .\\
\end{align}
In the final line, we're summing only over positive values of $k$, but we have two terms in the sum, one with a $+k$ and one with a $-k$.
Authors some times do this implicitly without telling you, which is confusing and stupid and on behalf of physicists everywhere I apologize for this sloppy terrible habit.
There is a nice reason to write things this way though.
We now have an integral that runs over each energy exactly once, and for each energy has two components: one move left and one moving right.
This may in some cases be a convenient way to write the general wavefunction, particularly if you know your particles are all coming from one direction and you want to write the wavefunction as a sum over energies.
Particle in a well
As you noted in the question, the approach to the bound particle already makes sense.
Fin

What I cannot seem to be able to understand is the seemingly different approach to the very same problem ... I just cannot reconcile the different approaches

It's just authors being sloppy.

$[a]$: For whatever it's worth, this exact issue confused me so much that I eventually stopped paying any attention to the details of any book's notation and just used them as a guide while developing my own notes and notation.
A: I will turn my comments into an answer:
The mathematics of solutions of Schrodinger equations were/are used to create models for specific data.
The k in the exponential describes momenta, it makes no sense in a plane wave, where energy is supposed to start at an x position and move to infinity, to have the sum  of two exponentials as representing probability amplitude , if we are trying to model a particle of energy E, a particle cannot have both a positive and a negative momentum, at least in time. I suppose that is (time) what drives the choice. Time cannot be ignored for real particle modelling.
In the infinite potential well solutions, the model is used to describe quantized states for the particle, where only E is a measurable variable to be modeled. There is no momentum that can be measured between bounds in quantum mechanics. That is why in the hydrogen atom we have orbitals and not orbits. So the general solution can be used, no need for a traveling plane wave.
I am attempting to say that it is the data which drive the choice of what solutions or combination of solutions to use to describe a specific physical set up. 
A: It some sense, it is unfortunate that the first two problems that many students encounter are the free particle and the infinite square well, because both problems are quite subtle from a technical standpoint.
Let's take a look at your questions first:

then they claim that k is a real non-negative number (why? I have no idea)

There's no reason $k$ should be non-negative - DanielSank mentions this thoroughly in his answer.  It must, however, be real, and this forces the energy to be non-negative.  

On the other hand, when the text go on to solve infinite potential well, they at first glance do something completely different, despite the problem, within the well, is the very same

The problems are not the same - the difference in domain radically alters the nature of the problem.  In the case of the infinite square well, the Hamiltonian has a discrete spectrum, and therefore a complete set of eigenvalue/eigenstate pairs, from which one can construct any arbitrary function in the Hilbert space.
On the other hand, in the case of the free particle, the Hamiltonian has a continuous spectrum.  This means that it does not have any eigenstates in the Hilbert space, so the above recipe fails spectacularly.  It does, on the other hand, yield a set of unphysical quasi-states, which do not live in the Hilbert space, but from which we can construct the physical solutions which do.  From there, we can look at the time-evolution of the physical states by examining the time-evolution of the unphysical ones from which it is built.
I will skim over each problem, and give you some context for why one might take slightly different approaches to them.  I will be somewhat technical, but hopefully not overly so.

When you write down a problem in quantum mechanics, you need to define a Hamiltonian operator and a Hilbert space of states on which it may act.  If the system corresponds to some subset of physical space (say, a particle on a line) then a good guess for the underlying Hilbert space is $L^2(\mathbb R)$, which essentially consists of all functions $f:\mathbb R \rightarrow \mathbb C$ such that
$$\int_{\mathbb R} |f(x)|^2 dx < \infty$$
equipped with the inner product
$$ \langle f,g\rangle  = \int_\mathbb R \bar f(x)\cdot g(x) \ dx$$
This is a good choice because we interpret $\int_a^b |f(x)|^2 dx$ as being proportional to the probability that the particle is measured to have a position in the interval $[a,b]$ (with equality holding if $f$ is normalized).  We should at the very least demand that a state be normalizable - otherwise this interpretation falls apart.
Once we have decided on the Hilbert space, we need a Hamiltonian operator $\hat H$ which is self-adjoint on our Hilbert space.  Notice that $\hat H$ needs to be self-adjoint, not merely Hermitian.  Those words are often taken to be synonyms, but they are not.  The distinction has to do with the domains on which unbounded operators such as $\hat X,\hat P,$ and $\hat H$ are allowed to act; I won't actually carry out any such calculations, but I'll mention when operators are Hermitian but not self-adjoint.  
Anyway, once we decide on a Hamiltonian, we typically seek to find its eigenstates, from which we can build our general solution. This is not always possible (e.g. in the case of a free particle), but it's a reasonable starting point.  Let's take a look at the actual problems in your question.

Free Particle
A free particle is just a particle on a line, so our Hilbert space should be $L^2(\mathbb R)$.  The momentum operator $\hat P$ operates on a state like this:
$$\hat P \psi := -i\hbar \psi'$$
and the Hamiltonian operator is simply
$$\hat H\psi := \frac{1}{2m}(\hat P \circ \hat P) \psi = -\frac{\hbar^2}{2m} \psi '' $$
We seek eigenfunctions of the Hamiltonian operator - that is, we seek functions $f\in L^2(\mathbb R)$ such that
$$ -\frac{\hbar^2}{2m} f''= E f $$
$$\implies f'' = -\frac{2mE}{\hbar^2} f$$
for some eigenvalue $E$.  There are three non-trivial subtleties here:


*

*An arbitrary element $f\in L^2(\mathbb R)$ is not even continuous, much less twice differentiable, and

*Even if $f$ is twice differentiable, there is no guarantee that $f''\in L^2(\mathbb R)$, and

*The two linearly independent candidates, namely the complex exponentials $e^{\pm ikx}$, are not even in $L^2(\mathbb R)$ themselves


The first two issues can be swept under the rug by only considering functions $f\in L^2(\mathbb R)$ which are twice differentiable, and whose second derivatives are also in $L^2(\mathbb R)$ (the set of such functions is called the Hilbert-Sobolev space $\mathcal H^2(\mathbb R)$).
The third issue is more problematic.  The typical way to handle it is to treat the complex exponentials as "quasi-states" which are unphysical, but still satisfy the above differential equation.  If $f=\exp(ikx)$, then
$$f'' = -k^2 f = -\frac{2mE}{\hbar^2}f$$
We then ask whether we could build a state in $L^2(\mathbb R)$ out of such states.  We can, and the result is a Fourier transform:
$$ f = \int_{\mathbb R} g(k) e^{-ikx} dk $$
From Parseval's theorem, 
$$\int_{\mathbb{R}} |f|^2 \ dx = \int_{\mathbb R} |g|^2 \ dk$$
So as long as $\int_{\mathbb R}|g|^2 \ dk < \infty$ then we are guaranteed that the resulting $f\in L^2(\mathbb R)$, and therefore corresponds to a physical state.  
The time dependent state is then
$$ f(x,t) = U(t,0) f(x,0) = e^{-i\hat H t} \int_{\mathbb R} g(k) e^{-ikx} dk$$
$$ = \int_{\mathbb R} g(k) e^{-ikx} e^{-i \frac{\hbar^2 k^2}{2m} t} \ dk$$
Notes:


*

*The plane waves $e^{-ikx}$ are not solutions to the time-independent Schrodinger equation, because the TISE is an eigenvalue equation on the Hilbert space $L^2(\mathbb R)$, and the complex exponentials do not belong to that space.

*Even so, we can build elements of $L^2(\mathbb R)$ from these unphysical quasi-states via superposition, which takes the form of a Fourier transform over momentum space; the resulting physical states evolve just like a superposition of the unphysical ones.

*The TISE does not have any eigenfunctions or eigenvalues in $L^2(\mathbb R)$, which means there is no such thing as physical states with definite energy for this system.


The last point is crucial and often slightly overlooked. We often claim that self-adjoint operators have complete sets of eigenstates, but this is generally not true if the spectrum of the operator is continuous, as it is here.
As an aside - it's easy to show that $\hat H$ is Hermitian, but somewhat tricky to show that it self-adjoint, which is a stronger requirement.  The same is true of $\hat P$ - it is self-adjoint (and therefore Hermitian), but it has a continuous spectrum, and so it has no physical eigenstates (though as before, we can treat the complex exponentials as quasi-states of definite momentum).

Infinite Square Well
In this problem, the particle is artificially restricted to lie in the interval $I=[0,L]$.  The Hilbert space $L^2(\mathbb R)$ is no longer an acceptable choice, so we make the modification that the Hilbert space underlying our system is now $L^2(I)$.
We actually need to make another restriction - we demand that physical states $\psi$ obey the boundary conditions $\psi(0)=\psi(L)=0$.  Therefore, our true Hilbert space $h$ (I'm sorry, we're running out of h's) is given by
$$h :=\{\psi \in L^2(I) \big| \psi(0)=\psi(L)=0\}$$
As before, we take the Hamiltonian operator to be 
$$ \hat H \psi := -\frac{\hbar^2}{2m} \psi ''$$
which we permit to act on functions which live in the domain $\mathcal D_H$, where
$$\mathcal D_H := \{\psi \in \mathcal H^2(I) \big| \psi(0)=\psi(L)=0\}$$
Unlike the previous case, the TISE actually yields solutions which belong to our Hilbert space, and which take the form
$$\psi_n = c_n \sin\big(\frac{n\pi x}{L}\big)$$
with eigenvalues
$$E_n = \frac{n^2\pi^2\hbar^2}{2mL^2}$$
Notes:


*

*As before, the Hamiltonian operator is self-adjoint on this Hilbert space.  However, since the spectrum of $\hat H$ is purely discrete, it does have a complete spanning set of eigenstates, as well as a notion of physical states of definite energy.

*One might imagine that the same is true of $\hat P$.  However, it can be shown that $\hat P$ is not self-adjoint on this Hilbert space.  Therefore, not only is there no such thing as a state of definite momentum, there isn't even a well-defined notion of what momentum is.  


Note that $\hat P$ as defined above actually takes elements of $h$ out of the Hilbert space.  This is due to our extremely restrictive boundary condition that $\psi(0)=\psi(L)=0$.  If we loosen it to periodic boundary conditions, so that $\psi(0)=\psi(L)$ but they are not necessarily zero, then it turns out that we can create a self-adjoint momentum operator.  This corresponds to a particle on a ring, and is an important problem to look at later.
