Skin grafting is a common technique employed to treat patients after burn injuries. In contrast to the frequency and gravity of contractures following from skin grafts, the phenomenon itself is still poorly understood and subject of studies. The development of an accurate model of skin contractures will allow medical researchers to better understand the healing process compared to current in vitro experiments and thus enable the design of efficient and patient-specific treatments. In this work we provide a mathematical model able to capture both the mechanical and biological processes involved in skin graft healing. Two-dimensional morphoelastic equations are used to model the mechanics of contraction of skin during the healing process. A system of equations describing the cell components (myo-, fibroblasts, collagen density and signalling molecule) is solved separately to derive the forcing terms for the mechanical system. The finite element method with bilinear quadrilateral elements is used to solve the differential equations on a moving time-dependent domain. A flux corrected transport (FCT) algorithm is used to stabilize the biological equations, thus enforcing positivity of the constituents in the dermal layer and nonlinearities are treated using Picard iterations. Time stepping is performed by applying the Euler backward scheme for the mechanical system and the implicit midpoint rule for the biological system. The results are replicated using the IGA (Isogeometric Analysis) library G+Smo in C++. The results of the IGA implementation provide a solid basis for future algorithms using IGA.