The statement is true if $χ_2^2+\frac{(χ_1χ_3)^2}{χ_1^2+χ_3^2} ≥ χ_4^2/2$, otherwise it might be false.
First, let $c = b-a$ so that the equation to solve becomes
$$\begin{align*}
xx' + yy' + zz' &= χ_1^2cc' + χ_2^2Λcc'Λ + χ_3^2(I_n-Λ)cc'(I_n-Λ)+χ_4^2(a(c+a)'Λ+Λ(c+a)a')\\
&= (χ_1^2+χ_3^2)cc' + (χ_2^2+χ_3^2)Λcc'Λ - χ_3^2(Λcc' + cc'Λ) + 2χ_4^2aa' + χ_4^2 (ac'Λ+Λca') \text{.}
\end{align*}$$
We will guess that $x,y,z$ are of the form $α_ic + β_iΛc + γ_ia$ for $i=1,2,3$, since this is the simplest way to build vectors out of the data we are given. Then
$$\begin{align*}
xx' + yy' + zz' &= \sum_{i=1}^3 α_i^2 cc' + β_i^2 Λcc'Λ + γ_i^2aa' + α_iβ_i(cc'Λ+Λcc') + α_iγ_i(ca'+ac') + β_iγ_i(Λca'+ac'Λ)\\
&= α^2 cc' + β^2 Λcc'Λ + γ^2aa' + αβ(cc'Λ + Λcc') + αγ(ca'+ac')+βγ(Λca'+ac'Λ) \text{.}
\end{align*}$$
(I write $αβ$ for the scalar product of $α$ and $β$.) Matching this expression against the target expression gives the following system.
$$\begin{align*}α^2 &= χ_1^2+χ_3^2\\
β^2 &= χ_2^2+χ_3^2\\
γ^2 &= 2χ_4^2\\
αβ &= -χ_3^2\\
αγ &= 0\\
βγ &= χ_4^2
\end{align*}$$
Let's suppose we are not in an edge case with $χ_i = 0$. We can solve this system geometrically. Choose two orthogonal vectors $α$ and $γ$ with the desired lengths. Then $αβ = -χ_3^2$ and $βγ = χ_4^2$ tell us the projection of $β$ on the plane generated by $α$ and $γ$. The minimal square-norm of $β$ is (I skip the computation) $χ_4^2/2 + \frac{χ_3^4}{χ_1^2+χ_3^2}$. Since $β$ should have a square-norm of $χ_2^2+χ_3^2$, in order to have a solution, we must have $χ_2^2+χ_3^2 ≥ χ_4^2/2 + \frac{χ_3^4}{χ_1^2+χ_3^2}$, or in other words $χ_2^2+\frac{(χ_1χ_3)^2}{χ_1^2+χ_3^2} ≥ χ_4^2/2$. If this is the case, there is a solution.
To find a counter-example, we need to take $χ_2^2+\frac{(χ_1χ_3)^2}{χ_1^2+χ_3^2} < χ_4^2/2$. For instance, with $n=3$ and $u=1$, take $c = (1,1,1)$, $a=(1,2,0)$, $χ_1 = 1$, $χ_2 = 0$, $χ_3 = 1$, $χ_4 = 3$. Then if my computations are correct, we find
$$\begin{pmatrix}
37&64&1\\
64&109&1\\
1&1&2
\end{pmatrix}$$
which is not positive semidefinite.
The condition $χ_2^2+\frac{(χ_1χ_3)^2}{χ_1^2+χ_3^2} ≥ χ_4^2/2$ might even be necessary. To show that, one could try to use the fact that even if this condition fails, there is a solution $β$ but with complex coefficients. This can be used to build a solution $x,y,z$, but there are complex coefficients: can this be used in some way to deduce that $xx'+yy'+zz'$ is not positive definite?
PS: Thank you for the context.
Let $(\psi,B)=(\chi',A')$; then \begin{gather*}
\psi'+2B\psi-\psi^2=12 \\
B'-B\psi+2B^2=24 \\
B^2-\psi^2=12 \tag{1}
\end{gather*} We can try to remove some nonlinearities by taking linear combinations. Adding, $$\psi'+B'+B\psi+2B^2-\psi^2=36$$ Now substituting the hyperbola condition (1): $$\psi'+B'+B(\psi+B)=24$$
But at this point, our equation is secretly 1-D. Note that $$2B=(B-\psi)+(B+\psi)=\frac{12}{B+\psi}+B+\psi\tag{2}$$ Thus if $u=B+\psi$, then we have the following autonomous ODE: $$24=u'+\frac{1}{2}\left(\frac{12}{u}+u\right)u=u'+6+\frac{u^2}{2}$$ Canceling, we discover a disguised Riccati equation $$36=2u'+u^2$$
Riccati equations have a trick: let $u=\frac{2v'}{v}$; then $$36=\frac{4(v''v-v'^2)+4v'^2}{v^2}=\frac{4v''}{v}$$ Thus $$v=\alpha\sinh{\sqrt{9}t}+\beta\cosh{\sqrt{9}t}$$
Correspondingly, $$\psi+B=u=\frac{2v'}{v}=6\cdot\frac{\alpha\cosh{3t}+\beta\sinh{3t}}{\alpha\sinh{3t}+\beta\cosh{3t}}$$ Solving for $B$ and $\psi$ via (2) and matching your boundary conditions is something I leave to you.
Best Answer
This is actually a linear first-order differential system with variable coefficients, in terms of $x_3 = x$. The present system can be recast as $ \boldsymbol{\chi}'(x) = {\bf A}(x)\, \boldsymbol{\chi}(x) $, where $\boldsymbol{\chi} = (\chi_1, \chi_2, \chi_3, \chi_4)^\top$ and $$ {\bf A} = \begin{pmatrix} \Delta_1 & 0 & \Delta_2 & 0\\ 0 & \Delta_1 & 0 & \Delta_2\\ \Delta_3 & 0 & \Delta_4 & 0\\ 0 & \Delta_3 & 0 & \Delta_4 \end{pmatrix} . $$ The solutions are given by $\boldsymbol{\chi}(x) = e^{{\bf B}(x)}\, \boldsymbol{\chi}_0$, where $\boldsymbol{\chi}_0 = (1,0,0,1)^\top$ and where $e^{{\bf B}(x)}$ denotes the matrix exponential of ${\bf B}(x) = \int_0^x {\bf A}(\xi)\, \text d \xi$, i.e. the matrix exponential of $$ {\bf B} = \begin{pmatrix} d_1 & 0 & d_2 & 0\\ 0 & d_1 & 0 & d_2\\ d_3 & 0 & d_4 & 0\\ 0 & d_3 & 0 & d_4 \end{pmatrix}, \qquad d_n(x) = \int_0^x \Delta_n(\xi)\,\text d \xi\, . $$ Using symbolic calculus software, the following $x$-dependent coefficients $\chi_n$ are obtained: $$ \boldsymbol{\chi} = \frac{e^{(d_1 + d_4)/2}}{2\sqrt{\mathcal{D}}} \begin{pmatrix} (d_1-d_4+\sqrt{\mathcal{D}})e^{\sqrt{\mathcal{D}}/2} - (d_1-d_4-\sqrt{\mathcal{D}})e^{-\sqrt{\mathcal{D}}/2} \\ 2\, (e^{\sqrt{\mathcal{D}}/2} - e^{-\sqrt{\mathcal{D}}/2})\, d_2 \\ 2\, (e^{\sqrt{\mathcal{D}}/2} - e^{-\sqrt{\mathcal{D}}/2})\, d_3 \\ (d_1-d_4+\sqrt{\mathcal{D}})e^{-\sqrt{\mathcal{D}}/2} - (d_1-d_4-\sqrt{\mathcal{D}})e^{\sqrt{\mathcal{D}}/2} \end{pmatrix} $$ where $\mathcal{D} = (d_1 - d_4)^2 + 4d_2d_3$. The corresponding Maple code is given below: