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. [1], 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) = ξ T 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 ∈ [−1, +1]. The procedure leads to a set of linear equations in terms of the unknown coefficients vector, ξ, 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.