Methods for fitting nonhomogeneous Markov models to panel observed data using direct numerical solution to the Kolmogorov Forward equations are developed. Nonhomogeneous Markov models occur most commonly when baseline transition intensities depend on calendar time, but may also occur with deterministic time-dependent covariates such as age. We propose transition intensities based on B-splines as a smooth alternative to piecewise constant intensities and also as a generalization of time transformation models. An expansion of the system of differential equations allows first derivatives of the likelihood to be obtained which can be used in a Fisher scoring algorithm for maximum likelihood estimation. The method is evaluated through a small simulation study and demonstrated on data relating to the development of cardiac allograft vasculopathy in post-transplantation patients.