For differential equations, the formal adjoint sets up an efficient way to integrate by parts, whether in one dimension or in many dimensions. For example,
when you integrate by parts, derivatives are stripped from the left and are moved to the right with a negative.
$$
\int \left(a\frac{d}{dx}f\right)gdx =-\int f\left(\frac{d}{dx}(ag)\right)dx
$$
So the "formal adjoint" of $a\frac{d}{dx}$ is $-a'+a\frac{d}{dx}$. This is a very old idea that goes back to Lagrange who then noted that the difference is an exact derivative:
$$
\left(a\frac{df}{dx}\right)g+f\left(\frac{d}{dx}(ag)\right)
= \frac{d}{dx}(afg).
$$
Lagrange used this to come up with reduction of order for ODEs, and other brilliant techniques that Lagrange was inclined to invent. You can do something similar with any ODE, and with some PDEs. For ODEs you end up with an integration identity
$$
\int (Lf)gdx - \int f(L^{\dagger}g)dx = \int \frac{d}{dx}(\cdots)dx = (\cdots)+C.
$$
If you carefully defined endpoint conditions on intervals, or boundary conditions for PDEs, you then get an operator relation
$$
(Lf,g) = (f,L^{\dagger}g).
$$
It's generally a lot of work and tedious detail to end up with a true operator adjoint in the strictest Functional Analysis sense. The adjoint in an inner product space is defined on the domain $\mathcal{D}(L^{\star})$ consisting of all $g$ for which there exists $h$ such that
$$
(Lf,g) = (f,h),\;\;\; \forall f\in\mathcal{D}(L).
$$
If $L$ has a dense domain, then you up with a unique $h$, which is defined as $L^{\star}g$, for which $(Lf,g)=(f,L^{\star}g)$.
As an example, consider $Lf=-f''$ on $C^2$ functions in $L^2[0,1]$. The formal adjoint of $L$ is $L^{\dagger} f=-f''$. The Functional Analysis adjoint is $L^{\star}f=-f''$ on the domain consisting of twice absolutely continuous on $[0,1]$ for which $f(0)=f'(0)=f(1)=f'(1)=0$. On this domain
$$
(Lf,g) = (f,L^{\star}g),\;\;\; f\in\mathcal{D}(L),\; g\in\mathcal{D}(L^{\star}).
$$
If you add a condition to the domain of $L$ (which has nothing to do with formal adjoint) such as $f(0)=f'(1)=0$, then $L^{\star}$ has domain consisting of twice absolutely continuous functions on $[0,1]$ for which the complementary conditions $f'(0)=f(1)=0$ hold. Again you get $(Lf,g)=(f,L^{\star}g)$ for $f\in\mathcal{D}(L)$ and $g\in\mathcal{D}(L^{\star})$. If you put all endpoint conditions on $\mathcal{D}(L)$, then the adjoint $L^{\star}$ consists of all twice absolutely continuous functions, but with no endpoint conditions.
There's a balancing act to get selfadjoint, where the conditions on $L$ must be chosen just right so that $L^{\star}$ has the some conditions on its domain. And the domain of $L$ must enlarged from twice continuously differentiable functions to twice absolutely continuous functions instead. When you then move to PDEs the domain of the adjoint will include all weak derivatives, which are no longer easily characterized in terms of absolute continuous. You find yourself in the realm of Sobolev spaces in order to have truly selfadjoint operators. You can see how the adjoint requirement $(Lf,g)=(f,h)$ for some $h$ and all $f\in\mathcal{D}(L)$ becomes a requirement of weak derivatives of some kind.
So, formal adjoint really just means that you're just manipulating differential operators and ignoring boundary conditions, or you're working in a classical setting of continuously differentiable functions, and haven't moved yet to a full Functional Analysis requirement of adjoint, closed operators, and abstract boundary conditions.
Laplacian Example: Suppose you have a nice bounded region $\Omega\subset \mathcal{R}^3$ with smooth boundary. The Laplacian $\nabla^2$ is formally selfadjoint because
$$
(\nabla^{2}f)g - f(\nabla^{2}g) = \nabla\cdot\{ g\nabla f - f\nabla g\}.
$$
Integrating over $\Omega$ ($dV$ is volume measure, $dS$ is surface measure) gives
$$
\int_{\Omega}(\nabla^2 f)g dV = \int_{\Omega}f(\nabla^2 g)dV
= \int_{\partial\Omega} g \frac{\partial f}{\partial n}-f\frac{\partial g}{\partial n} dS.
$$
If you let $f \in \mathcal{C}_c^{\infty}(\Omega)$, then $f$ and $\frac{\partial f}{\partial n}$ vanish on $\partial\Omega$, which means that $g$ can be unrestricted on $\partial\Omega$. So that's not going to lead to an actual selfadjoint problem. That's the same as looking at $\frac{d^2}{dx^2}$ on $[0,1]$ and assigning $f(0)=f'(0)=f(1)=f'(1)=0$, which is the case I mentioned above. However, if the functions $f,g$ are $C^2$ right up to the boundary and vanish on $\partial\Omega$, then that's a more symmetric condition, and you can see that $(\nabla^2f,g)=(f,\nabla^2g)$. An alternative condition is where $\frac{\partial f}{\partial n}=\frac{\partial g}{\partial n}=0$ vanish on the boundary, which is also a symmetric type of condition.