This study is devoted to designing two hybrid computational algorithms to find approximate solutions for a class of singularly perturbed parabolic convection–diffusion–reaction problems with two small parameters. In our approaches, the time discretization is first performed by the well-known Rothe method and Taylor series procedures, which reduce the underlying model problem into a sequence of boundary value problems (BVPs). Hence, a matrix collocation technique based on novel shifted Delannoy functions (SDFs) is employed to solve each BVP at each time step. We show that our proposed hybrid approximate techniques are uniformly convergent in order
O
(
Δ
τ
s
+
M
−
1
2
)
{\mathcal{O}}\left(\Delta {\tau }^{s}+{M}^{-\tfrac{1}{2}})
for
s
=
1
,
2
s=1,2
, where
Δ
τ
\Delta \tau
is the time step and
M
M
is the number of SDFs used in the approximation. Numerical simulations are performed to clarify the good alignment between numerical and theoretical findings. The computational results are more accurate as compared with those of existing numerical values in the literature.