  In this paper a technique is suggested to integrate linear initial boundary
value problems with exponential quadrature rules in such a way that the order
in time is as high as possible. A thorough error analysis is given for both the
classical approach of integrating the problem firstly in space and then in time
and of doing it in the reverse order in a suitable manner. Time-dependent
boundary conditions are considered with both approaches and full discretization
formulas are given to implement the methods once the quadrature nodes have been
chosen for the time integration and a particular (although very general) scheme
is selected for the space discretization. Numerical experiments are shown which
corroborate that, for example, with the suggested technique, order $2s$ is
obtained when choosing the $s$ nodes of Gaussian quadrature rule.
