Can I state that $\Psi (x_1, \dots , x_n,t)= \sum_{i=1}^n a_i \psi (x_i,t) $ via superposition? Given that the hamiltonian $\hat H$ of a system is a linear operator and $\dot \psi (x_i,t)$ does not depend on spatial coordinates $x_1, ..., x_n$ with bases $\hat e_1, ... , \hat e_n$ can I state that $$\Psi (x_1, \dots , x_n,t)= \sum_{i=1}^n a_i \psi (x_i,t) $$ and $\sum_{i=1}^n |a_i|^2 = 1 $; which implies the superposition principle? I am developing an algorithm for solving Schroedinger Equations and I want to know if my approach, namely solving for each spatial dimension and adding a time-dependent phase, and then combining solutions via the superposition principle is a good approach to the problem.
If yes, are there other conditions that must be satisfied? 
$\psi (x_i, t) $ are the wavefunctions corresponding to each spatial dimension (which are as well time dependent) and the system contains one particle in $n$-dimensional space
 A: No, you can't do that.
First of all because it doesn't make much sense to talk about the wave function corresponding to a spatial dimension.
A wave function $\psi(x_1,\dots,x_n,t)$ gives you the probability amplitude of the system being at the position $(x_1,\dots,x_n)$ at the time $t$. It must by definition be dependent on all the coordinates necessary to describe the system (which can be spatial coordinates, spin degrees of freedom or whatever).
Even assuming that with $\psi(x_i,t)$ you mean some peculiar wave function which does depend on all the spatial coordinates $(x_1,\dots,x_n)$, but which is constant with respect to variations of $x_j, j \neq i$, this just doesn't work.
To see this consider the simple example of a free particle in the usual three dimensions. The plane wave solution is, neglecting normalization factors,
$$ \tag 1 \Psi(x_1,x_2,x_3,t) = e^{i(k_1 x_1+k_2 x_2+k_3 x_3)} e^{-i E_k t/\hbar}, $$
where
$$ E_k \equiv \frac{\hbar^2 k^2}{2m} $$
is the energy of the particle with wave vector $\mathbf k$.
You can easily see that (1) cannot be decomposed in a sum of functions depending only on a single spatial variable, indipendently from the time dependence factor.
If however you were asking about decomposing the wave function as a product of wave functions depending on the single coordinates, than the answer would be almost.
It is a general feature that when two degrees of freedom are not interacting, the correspoding time-independent wave function can be factored as a product of time-independent wave functions each depending on one degree of freedom.
The time-dependence is applied to the whole product though, not to the single factors.
An example of this is again (1), as you can easily verify.
This factorization can be done when the Hamiltonian is a sum of terms each acting on a single coordinate. Something like, in your example,
$$ \tag 2 \mathcal H = \sum_{i=1}^n \mathcal H_i$$
where $\mathcal H$ is the total Hamiltonian and $\mathcal H_i$ are each acting on only $x_i$. Then you can solve each time-independent, one-dimensional Schrödinger equation, i.e. solve for the whole set of eigenfunctions $\psi_i^{(m)}(x_i)$ of the one-dimensional problem with Hamiltonian $\mathcal H_i$:
$$ \tag 3 \mathcal H_i \psi_i^{(m)}(x_i) = E_i^{(m)} \psi_i^{(m)}(x_i), $$
where $m$ labels the various solutions of (3).
Once you've done this for each $i$, you have a base of eigenfunctions in which you can decompose the total wave function $\Psi$. Namely, including also the time-dependence, the solution of the time-dependent Schödinger equation of the whole system can be decomposed as
$$ \tag 4 \Psi(x_1,\dots,x_n,t) = \sum_{\{ m_i \}_{i=1,\dots,n}} c(\{m_i \},t) \prod_{k=1}^n \psi_k^{(m_k)}(x_k). $$
Given any initial condition, your problem is reduced to finding the time-dependence of the coefficients $c(\{m_i\}_i,t)$, which represent the probabiliy amplitude of the system being in the state labeled by the set of quantum numbers $\{m_i\}_i$ at the time $t$.
As you can imagine, this is in most cases highly non trivial.
