This thesis presents the formulation and implementation details of a method based on the modified error in constitutive equation to address the viscoelasticity imaging problem in situations in which information on boundary conditions is unknown. In the viscoelasticity imaging problem, the goal is to produce images of the viscoelastic properties of a material from displacements measured in its interior. These measured displacements are always corrupted with noise, which poses a challenge to methods designed to solve such problems. Moreover, in practical applications such as biomedical imaging, where the material of interest is tissue, the magnitude and spatial distribution of the excitation used to generate the displacements that are measured are not exactly known. This is a challenge to optimization-based methods for the imaging problem, as the lack of boundary conditions leads to ill-posed forward problems. The method developed here overcomes this challenge and at the same time handles noisy and incomplete displacement data. This thesis is divided into two chapters. The first chapter, which is an adaptation of a journal article that has been submitted for publication, presents the method and relevant derivations. Results from numerical experiments are also included in this chapter. The second chapter details the implementation of the method in the DinamicaE simulation suite, developed by the Computational Mechanics and Inverse Problems Group led by Professor Wilkins Aquino. The DinamicaE simulation suite is the result of over four years of development effort, which is still ongoing. Development of a large component of this software has been one of the author's main contributions as a Ph.D. student. DinamicaE is a massively parallel research code that solves problems in steady-state dynamics, acoustics, and acoustic-structure interaction. Moreover, the software also solves imaging problems in these domains. One of the main goals of DinamicaE is to assess the feasibility of algorithms such as the one presented in this thesis to solve problems of interest in the field of biomedical imaging, which seeks to provide early diagnosis for many physical illnesses. Its modular design, which is made possible by the features offered by the C++ programming language, allows for simple implementations of the algorithm presented in this thesis and many more.