Finding the decomposition of a function on a cubic spline basis of functions In a computational project, I need to solve a partial differential equation. Standard procedure is to consider the weak formulation of the problem which maps it onto an algebraic problem. 
With cubic hermite splines,if my space is discretized into $n$ points, I can form a basis of $2n$ functions $\{\phi_i(x)\}$ to express any solution as a linear combination of these basis functions. Then, whatever function $f(x)$ is approximated by its projection onto that basis, that is to say
$$ f(x) \simeq proj_f(x) = \sum\limits_{i=1}^{2n} a_i \phi_i(x) $$
In order to verify whether my simulations work, I need to plug in an analytical solution g(x) into the machinery in terms of its decomposition on the basis. So, given $g(x)$ and $\{\phi_i(x)\}$, I want to find ($\vec{a}$), such that 
$$g(x) \simeq proj_g(x) \equiv \sum \limits_{i=1}^{2n}a_i\phi_i(x)$$
Is the following reasoning correct (?):
$$ \underbrace{\int \limits_x g(x)\phi_j(x) dx}_{b_j} \simeq \sum \limits_{i}^{2n}a_i\underbrace{\int \limits_x \phi_i(x)\phi_j(x) dx}_{M_{ij}} \quad \quad \forall j $$
So essentially, this gives me one equation for each j, so it is equivalent to 
$$ \vec{b} = M\cdot\vec{a} \quad \to \quad \vec{a} = M^{-1}\cdot\vec{b}$$
Now when I do this computationally, and retrieve $\vec{a}$, and plot the function $proj_g(x) \equiv \sum \limits_{i} a_i\phi_i(x)$, it doesn't look like the function $g(x)$. 
SO here is my question : is the reasoning above wrong ? or is it right and I need to look for an implementation error ?

The reasoning and formulas are correct in principle: you look for the least squares approximation to $g$, the orthogonal projection in the given inner-product space on the subspace generated by $\{\varphi_i\}$.
