Rigorous justification for this formal solution to $f(x+1)+f(x)=g(x)$ Let $g\in C(\Bbb R)$ be given, we want to find a solution $f\in C(\Bbb R)$ of the equation 

$$
f(x+1) + f(x) = g(x).
$$

We may rewrite the equation using the right-shift operator $(Tf)(x) = f(x+1)$ as 
$$
(I+ T)f=g.
$$ 
Formally, I can say that the solution of this equation is 

$$
f= (I+ T)^{-1} g.
$$

Of course, I am aware that there are infinitely many solutions to the equation since the kernel of $(I+T)$ consists of all $h\in C(\Bbb R)$ such that $h(x+1)=-h(x)$, e.g. $\sin(\pi x)$, but please bear with me for a moment here.
By the theory of operator algebra, if $f,g$ are from some nice Banach space $X$ and our linear operator $T:X\to X$ satisfies $\|T\|<1$, then we have 
$$
f = \left(\sum_{n=0}^\infty (-T)^n \right)g.
$$
However, it is not unreasonable to expect that we should have $\|T\|=1$ for a right-shift operator in most reasonable function spaces so let's try to solve the equation 
$$
f= (I+\lambda T)^{-1} g.
$$
for $|\lambda| <1$ first then we'll take $\lambda\to 1$. Note that all the steps until now is purely formal since $C(\Bbb R)$ is not a normed space.

To illustrate what I meant, let's say we take $g(x) = (x+2)^2$. We now try to implement the above method (for $|\lambda|<1$) to get
$$\begin{align}
f(x) &= \left(I - \lambda T + \lambda^2 T^2 - \dots \right) g(x) \\
&= (x+2)^2 -\lambda (x+3)^2 + \lambda^2 (x+4)^2 + \dots \\
&= \left(1-\lambda+\lambda^2-\dots \right)x^2 + \left(2-3\lambda+4\lambda^2-\dots \right)2x + \left(2^2-3^2\lambda+4^2\lambda^2-\dots \right) \\
&=  \frac{1}{1+\lambda} x^2 + 2 \frac{2+\lambda}{(1+\lambda)^2} x +  \frac{4+3\lambda + \lambda^2}{(1+\lambda)^3}.
\end{align}$$
We shall be brave here and substitute $\lambda=1$ even though the series doesn't converge there. This gives
$$
f(x) = \frac 12 x^2 + \frac 32 x + 1
$$
but voilà, for some mysterious reasons unknown to me, this $f$ actually solves our original equation $f(x+1) + f(x) = (x+2)^2$ !

My question is simply:

What are the hidden theories behind the miracle we observe here? How can we justify all these seemingly unjustifiable steps?

I can't give you a reference to this method because I just conjured it up, thinking that it wouldn't work. To my greatest surprise, the answer actually makes sense. I am sure that similar method is probably practiced somewhere, probably by physicists.
Remark: I posted an isomorphic question on MSE earlier and one of the commenters mentioned that this could be related to the holomorphic functional calculus (HFC), specifically the part where I let $\lambda \to 1$. I once learned HFC just for an exam and don't have much memory of it so I don't immediately see if we can make the above method fully rigorous using merely the standard HFC or not.
 A: If $g(x)$ grows slower than any exponential as $x$ goes to $\infty$, then $ \sum_n (-\lambda)^n T^n g$ clearly converges (absolutely) pointwise for any $\lambda <1$, and when $(I + \lambda T)$ is applied, recovers $g$.  
If this sequence happens to converge pointwise as $\lambda$ goes to $1$, then it will recover $g$ when $I+ T$ is applied. I don't know a great perfectly general criterion for this to converge pointwise, but it will for instance when $g$ is polynomial (or a sum of polynomials times exponential phases $e^{ \pi i t x}$ where $t$ is not an odd integer.)
In general, when there is not a unique solution, you can't expect these formal inverses to give a unique solution. But if anything that's even more reason to expect them to give you at least one solution.
From an analytic perspective, the proof that a formal solutions solves the equation by a telescoping sum is much simpler than the proof that it is the unique solution, and we should expect it to hold in more general settings. 
A: The "miracle" is just that $T$ is also defined $\mathbb{R}_n[X]\rightarrow \mathbb{R}_n[X]$, a elementary linear operator on finite dimensional space of polynomes. And on this set we have 
$$T=\begin{pmatrix} 1 & 1 & 1 & & \\ 0 & 1 & 2 & &\vdots \\ 0 & 0 & 1 & & \\ 0 & 0 & 0 & \ddots  \\  0 & 0 & 0 & 0 & 1 \end{pmatrix} $$
with $T_{ij}=\begin{pmatrix}i \\ j\end{pmatrix}$ if $i\geq j$ and $0$ otherwise. So $\det(I+\lambda T)=(1+\lambda)^{n+1}$ and then $(I+\lambda T)^{-1}$ is analytic on $\lambda\in \mathbb{C}/\{-1\}$.
A: A direct and precise way to arrive at your answer is to appeal to the theory of Fourier transforms of distributions. If I define the Fourier transform as
$$F(k)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty f(x)e^{ikx}\,dx,$$
then the equation $f(x+1)+f(x)=g(x)$ becomes upon Fourier transformation
$$e^{-ik}F(k)+F(k)=G(K)\Rightarrow F(k)=\frac{G(k)}{1+e^{-ik}}.$$
Inverse Fourier transformation gives the desired $f(x)$,
$$f(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty \frac{G(k)}{1+e^{-ik}}e^{-ikx}\,dk.$$
The Fourier transform of $g(x)=(x+2)^2$ is given in the distributional sense by derivatives of the Dirac delta function $\delta(k)$,
$$(2\pi)^{-1/2}G(k)=4  \delta (k)-4 i \delta '(k)-\delta ''(k).$$
These distributions are well defined upon integration, giving the inverse Fourier transform $f(x)=\frac{1}{2} \left(x^2+3 x+2\right)$.    
More general cases follow readily:
$$g(x)=(x+a)^2\Rightarrow f(x)=\tfrac{1}{2} (a+x-1) (a+x),$$
$$g(x)=(x+a)^3\Rightarrow f(x)=\tfrac{1}{8} (2 a+2 x-1) \left(2 a^2+a (4 x-2)+2 (x-1) x-1\right),$$
$$g(x)=(x+a)^4\Rightarrow f(x)=\tfrac{1}{2} (a+x-1) (a+x) \left(a^2+a (2 x-1)+(x-1) x-1\right).$$

A: Here is another approach, which like yours is not entirely rigorous, but I think gives a bit more insight into the situation.
If you are willing to assume that $f$ is not just smooth but also analytic, then by writing the $f(x+1)$ term as a Taylor expansion you can write the equation in the form
$$
(e^D+I) f = g
$$
where $D$ is the differentiation operator. It is then natural to apply the (formal) inverse operator and get that
$$
f = (e^D+I)^{-1} g.
$$
Now note that the function $t\mapsto (e^t+1)^{-1}$ itself has the Taylor expansion
$$
\frac{1}{e^t+1}
=\frac{1}{2}-\frac{t}{4}+\frac{t^3}{48}-\frac{t^5}{480}+\frac{17 t^7}{80640}+...
$$
This suggests that the solution to our equation can be expressed as
$$
f = \left(\frac{1}{2}I - \frac{1}{4}D + \frac{1}{48}D^3 - \frac{1}{480}D^5 +... \right)g.
$$
The nice thing about this expansion is that if $g$ is a polynomial then this is in fact a finite sum, so you get an expression for $f$ as a finite linear combination of $g$ and its derivatives. In your example of $g(x)=(x+2)^2$ we have
$$
f(x) = \frac12 g(x) - \frac14 g'(x) + \frac{1}{48} g'''(x)
$$
which indeed comes out to $f(x)=\frac12 x^2+\frac32 x + 1$.
So what is going on here? The situation is very analogous to the formal derivation of the Euler-Maclaurin summation formula, in which the reason for the appearance of the Bernoulli numbers in the formula can be explained by the fact that one is essentially inverting the operator $e^D-I$ (instead of $e^D+I$ in the current discussion). See this comment of Terry Tao's from his blog. Note that in that case as well the derivation is purely formal and one needs to do more careful work to verify that the resulting formula makes sense. 
Ultimately, I think what all of this illustrates is the principle that formal manipulations using operators involving infinite series, Taylor expansions etc (I think this goes under the name of operator calculus, loosely speaking, but perhaps someone else can give more precise terminology) often lead to correct results without there necessarily being some unifying principle at work that explains why this is happening (which is not to say that one shouldn't look for good explanations of course). I can't think of specific examples, but I think one runs into this phenomenon in lots of places in analysis, mathematical physics, probability and other places.
A: Just to add another answer to the pile...
I suspect that the 'magic' observed here may be related to so-called local spectral theory. It is possible to do quite a lot (at least in the Banach space setting) with the notion of a local resolvent $\rho_T(x)$ of a bounded Banach space operator $T:X\to X$; for $x\in X$, this is defined as the union of all open $U\subset\mathbb{C}$ for which there is some $f\in \mathscr{O}(U,X)$ such that $(T-\lambda)f(\lambda)=x$ for all $\lambda \in U$. The standard reference for this sort of stuff (at least for the Banach space setting) is the book of Neumann and Laursen. I don't know of a corresponding reference, but it wouldn't surprise me if versions of local spectral theory have been studied in the Frechet/barrelled settings as well.
