Given a linear second-order differential operator $${\mathcal {L}}\equiv \phi \,D^2+\psi \,D$$
L
≡
ϕ
D
2
+
ψ
D
with non zero polynomial coefficients of degree at most 2, a sequence of real numbers $$\lambda _n$$
λ
n
, $$n\geqslant 0$$
n
⩾
0
, and a Sobolev bilinear form $$\begin{aligned} {\mathcal {B}}(p,q)\,=\,\sum _{k=0}^N\left\langle {{\mathbf {u}}_k,\,p^{(k)}\,q^{(k)}}\right\rangle , \quad N\geqslant 0, \end{aligned}$$
B
(
p
,
q
)
=
∑
k
=
0
N
u
k
,
p
(
k
)
q
(
k
)
,
N
⩾
0
,
where $${\mathbf {u}}_k$$
u
k
, $$0\leqslant k \leqslant N$$
0
⩽
k
⩽
N
, are linear functionals defined on polynomials, we study the orthogonality of the polynomial solutions of the differential equation $${\mathcal {L}}[y]=\lambda _n\,y$$
L
[
y
]
=
λ
n
y
with respect to $${\mathcal {B}}$$
B
. We show that such polynomials are orthogonal with respect to $${\mathcal {B}}$$
B
if the Pearson equations $$D(\phi \,{\mathbf {u}}_k)=(\psi +k\,\phi ')\,{\mathbf {u}}_k$$
D
(
ϕ
u
k
)
=
(
ψ
+
k
ϕ
′
)
u
k
, $$0\leqslant k \leqslant N$$
0
⩽
k
⩽
N
, are satisfied by the linear functionals in the bilinear form. Moreover, we use our results as a general method to deduce the Sobolev orthogonality for polynomial solutions of differential equations associated with classical orthogonal polynomials with negative integer parameters.