  We present efficient MATLAB implementations of the lowest-order primal hybrid
finite element method (FEM) for linear second-order elliptic and parabolic
problems with mixed boundary conditions in two spatial dimensions. We employ
the Crank-Nicolson finite difference scheme for the complete discrete setup of
the parabolic problem. All the codes presented are fully vectorized using
matrix-wise array operations. Numerical experiments are conducted to show the
performance of the software.
