  This study shows how to obtain least-squares solutions to initial and
boundary value problems to nonhomogeneous linear differential equations with
nonconstant coefficients of any order. However, without loss of generality, the
approach has been applied to second order differential equations. The proposed
method has two steps. The first step consists of writing a constrained
expression, introduced in Ref. \cite{Mortari}, that has embedded the
differential equation constraints. These expressions are given in term of a new
unknown function, $g (t)$, and they satisfy the constraints, no matter what $g
(t)$ is. The second step consists of expressing $g (t)$ as a linear combination
of $m$ independent known basis functions, $g (t) = \mathbf{\xi}^T \mathbf{h}
(t)$. Specifically, Chebyshev orthogonal polynomials of the first kind are
adopted for the basis functions. This choice requires rewriting the
differential equation and the constraints in term of a new independent
variable, $x\in[-1, +1]$. The procedure leads to a set of linear equations in
terms of the unknown coefficients vector, $\mathbf{\xi},$ that is then computed
by least-squares. Numerical examples are provided to quantify the solutions
accuracy for initial and boundary values problems as well as for a control-type
problem, where the state is defined in one point and the costate in another
point.
