You will run into some issues with differential equations with singularities.
Consider the differential operator $x\frac{d}{dx} -t $ from the trivial rank $1$ vector bundle on $\mathbb R$ to itself, for some constant $t$. The adjoint map is $-xd/dx -1 - t$, which is another operator in the same class.
If $t$ is not a nonnegative integer, then this complex is locally exact. We have to check that the differential equation $xdy/dx - ty =f(x)$ has solutions for smooth $f$. We will check this for $t<0$, but I think it is also true when $t$ is not a nonnegative integer. A solution is:
$$y = \frac{ f(0) }{t} + x^t \int_{0}^x \left(f(z)-f(0)\right) z^{-1-t} dz $$
This gives a global solution also, so it is globally exact. Moreover, since any solution to the differential equation is a multiple of $x^{t}$, if $t$ is not a nonnegative integer, there are no nonzero solutions, so the regular cohomology is trivial.
If $t$ is a negative integer, then the dual complex will have nontrivial $H^0$ due to the nonzero solutions, otherwise, say for $-1<t<0$, we can easily show that the dual complex, with $t$ also in that range, will have nontrivial compactly supported $H^1$. Indeed, since there are no solutions to the homogeneous version of the differential equation, our solution of the inhomogeneous version is unique, and we can easily find a compactly supported $f(x)$ where the unique solution $y$ is not compactly supported.
On the other hand, suppose we have a smoothness condition - specifically, that the kernel of the first map is a locally constant sheaf. In other words, a local solution to the differential equations that define the first map can be extended uniquely along any path, with possible monodromy.
Given a differential equation, a common trick is to add enough extra functions to make the equation first order. We can just as easily do this with a complex of vector bundles with differential equation operators - add variables to each map in the reverse order. This process is a homotopy equivalence of complexes, as is its dual.
Take the locally constant kernel sheaf, view it as a vector bundle iwth flat connection, and tensor it with the de Rham complex. We will build a map from this complex to the original one. This is plausible because they are both injective resolutions of the same thing, but we need to check it can be done with vector bundle maps. This is trivial in degree $0$. If we have built maps for the first $n$ degrees, we compose the $n$th map with $d$ and get a first-order function on $\Omega^{n-1}$ that vanishes, locally, on the image of anything from $\Omega^{n-2}$. Such a function, by the linear algebra of differential forms at a single point, depends only on $d$ of the form on $\Omega_{n}$.
This bundle map is a quasi-isomorphism of sheaf complexes. If we can check that its dual is also a quasi-isomorphism, we win - duality in an arbitrary locally free complex can be reduced to duality for the de Rham complex. By using mapping cones, it is sufficient to check that if a first-order complex is locally exact, its dual is also locally exact.
Let $V_0 \to V_1 \to \dots V_n$ be a locally exact first-order complex. We will actually be able to find a homotopy to $0$. $d: V_0 \to V_1$ is a first-order differential equation with no local solutions. If it has no solutions,it must have a formal reason. Specifically, if $f_1,..f_k$ are local coordinates for $V_0$, then by taking linear combinations of the differential equations, their derivatives, and the commutation relations,we must be able to obtain $f_1,\dots,f_k$. Otherwise we could solve it along curves and extend consistently to the whole space.
But this linear combination just gives an operator $k: V_1 \to V_0$ such that $kd$ is the identity. Now we have a differential equation $k\oplus d$ on $V_1$, still linear, that has no solutions. Repeating the process, we eventually get a homotopy between the bundle and $0$. Applying this homotopy to the dual, it will be locally exact as well.
Best Answer
I'll have more time to write and provide a more thorough answer later, but I think the most straightforward proof (which I agree is hard to find) comes via sheaf theory: On the one hand, there is a sheaf of locally finite singular chains whose hypercohomology is your $H^{lf}$. I work out the details in the setting of intersection homology on pseudomanifolds in Section 3 of the following paper, but the ordinary (non-intersection) lf homology on locally compact (paracompact?) spaces is essentially a special case via basically the same arguments: "Singular chain intersection homology for traditional and super-perversities" Transactions of the American Mathematical Society 359 (2007), 1977-2019
On the other hand, you can find a good treatment of the dualizing complex in Borel's book on intersection homology, section V.7. The idea then is to show that these two sheaves are quasi-isomorphic. This is the part I'll have to think about and add latter, though I'll note that it's not so hard for manifolds: in that case, it's not so hard to show by hand that both have to be quasi-isomorphic to the orientation sheaf (shifted appropriately). That's shown in Borel for the dualizing complex, and follows for the sheaf of locally finite singular chains via the local homology computation on the manifold.
Sorry I don't have time to write more now, but hopefully that's a start on the ideas.
Added later:
Oookay, this was harder than I thought. I’m not sure I can give you a really satisfying answer, but I think I can convince you that what you want is contained in Bredon as Corollary V.12.21. This corollary pretty much says that if $X$ is nice enough (though still pretty general), then ${}_SH_p^{\Phi}(X;\mathcal B)\cong H_p^{\Phi}(X;\mathcal B)$, where here $\mathcal B$ is a sheaf and $\Phi$ is a paracompactifying family of supports. For what we want, let’s take $\Phi$ to be all the closed sets, so we need to assume $X$ paracompact. We’re also going to let $\mathcal B$ be the constant sheaf with stalk the ground ring $L$. Let’s call that $\mathcal R$ (since $\mathcal L$ is used for lots of other things in these sections of Bredon). In all of the definitions we’ll talk about, $\mathcal B$ would come in as a tensor product factor, so letting it be constant basically lets us ignore those factors.
Now, on the left we have ${}_SH_p^{cld}(X;\mathcal R)$, which is Bredon’s singular homology with (now constant) sheaf coefficients. He gives a pretty good explanation for this back on pages 287-289. The key point connecting this to what I wrote previously is the following. Suppose you start with the (homologically graded) presheaf he calls $\Delta_*^c(X,X-\bar U;R)$, which is just the usual relative singular chain complex. Since our spaces are locally compact, this sheafifies to the same sheaf as $\Delta_*^{lf}(X,X-\bar U;R)$, the sheaf of locally finite singular chains, whose hypercohomology (or hyperhomology with this indexing) computes the locally finite singular homology (more details about this can be found elsewhere in Bredon or in the Transactions paper of mine that I already mentioned). Though actually here you don’t even need the hypercohomology; on the top of page 289 he just takes the homology of the complex of sections (which is the same as the hyperhomology because the sheaf is homotopically fine – shown earlier in the book). The argument on page 288-289 further shows that the resulting homology groups are isomorphic to singular homology as he defined it on page 288, using basically the same sheaf but also employing the direct limit over subdivisions. So, the upshot is that the term on the left is locally finite singular homology.
For the term on the right, we need to trace through the definitions starting on page 292. By the boxes at the top, $$H_p^{cld}(X;\mathcal R)= H_p(C_*^{cld}(X;\mathcal R))= H_p(\Gamma(\mathcal C_*(X;\mathcal R))) =H_p(\Gamma(\mathcal D_*(\mathcal J^*(X;L)))\otimes \mathcal R),$$ where $L$ is our ground ring. So if $X$ is locally connected, we can ignore the $\mathcal R$ and we’re just taking the homology of the global sections of $\mathcal D_*(\mathcal J^*(X;L))$. I claim this is just the Verdier dualizing sheaf with homological indexing. To see this, we unwrap more definitions! By the top of page 292, $\mathcal J^*(X,L)$ is just the canonical injective resolution of the constant sheaf (he also calls the constant sheaf $L$ and I have called it $\mathcal R)$. Let’s just call this injective (and hence also c-soft) resolution $\mathcal K^*$. So we need to understand $\mathcal D(\mathcal K^*)$. By the definitions on pages 289-290, $\mathcal D$ should take a second argument $M$, but there’s a note on page 290 saying it’s omitted from the notation if $M=L$. So we’re really after $\mathcal D(\mathcal K^*;L)$.
By the second box on page 290, $\mathcal D(\mathcal K^*;L)$ means $\mathcal D(\Gamma_c\mathcal K^*;L)$, where $\Gamma_c\mathcal K^*$ is the cosheaf given by $(\Gamma_c\mathcal K^*)(U)=\Gamma_c(\mathcal K^*|U)$ by the definition part of Proposition V.1.6. But now piecing together the meaning of Definition V.2.1, which is given explicitly on the top of page 291, $\mathcal D(\Gamma_c\mathcal K^*;L)(U)=Hom(\Gamma_c\mathcal K^*(U),L^*)$, where $L^*$ is an injective resolution of $L$ according to the beginning of Section V.2. So, altogether, $\mathcal D(\Gamma_c\mathcal K^*;L)$ is just the sheaf $U\to Hom(\Gamma_c(\mathcal K^*|U),L^*)$.
But this is exactly the definition of the dualizing sheaf given by Borel on page 119 of Intersection Cohomology. That this is $f^!$ of the constant sheaf on a point is on Borel page 131.
So there you go. Tracing through the Bredon definitions is a nightmare, and not a satisfying conceptual thing at all, but it’s a proof.
To say just a couple of words about the proof of the Corollary, the point is that he defines at some point the idea of a ``flabby quasi-n-coresolution of L,’’ where here $L$ is the notation for the constant sheaf again. The point of these is that if you have one, say $\mathcal L_*$, and want to compute $H_*(X;\mathcal B)$ for some sheaf $\mathcal B$, then instead of going through the horrible definitions we just unwound, you have $H_*(X;\mathcal B)$ is just $H_*(\Gamma(\mathcal L_*\otimes \mathcal B))$. That’s Theorem V.12.20, and this is sort of like what happens in Borel V.6.7 (or other places) when you make a sheaf complex c-soft by tensoring it with a flat c-soft resolution of the constant sheaf. But a bit earlier in Theorem V.12.14 it’s shown that the same sheaf used to compute singular homology (the one using the subdivisions, called $\mathfrak S_*$) is such a flabby quasi-n-coresolution of $L$. That actually just follows from the definitions and the properties on the space. So it follows that $H_*(X;\mathcal B)\cong H_*(\Gamma(\mathfrak S_*) \otimes \mathcal B))$, which is the definition of singular homology with coefficients in $\mathcal B$.
It’s really a shame this isn’t all brought out more clearly in Bredon’s exposition.
Whew.