This contribution considers the use of Mathematica (Wolfram Research Inc.) to estimate the parameters of linear mixed-effects (LME) models by maximum likelihood (ML) and restricted likelihood (REML) methodologies. Mixed-effects models are currently applied to estimate population pharmacokinetic parameters from clinical data. Contrary to their fixed-effects counterparts, which simply address within-patient variability, these models take into consideration inter- as well as intra-individual variability. The model parameters are estimated by maximizing an approximation to the likelihood function using techniques such as Lindstrom and Bates algorithm, Laplacian and adaptive Gaussian methods. Considerable efforts have been devoted to integrating these procedures into computer programs (e.g., the NLMIXED procedure in SAS, LME and NLME in S-PLUS and R) to reach a broader scientific audience and facilitate the incorporation of these tools in routine clinical pharmacokinetic data analysis. The codes, offered in this work, were tested against two benchmark problems. The results were similar to those generated by other programs.