This paper presents a new algorithm for computing stable discrete minimal surfaces bounded by fixed or free boundary curves in 3D space, as well as in the 3-dimensional sphere and hyperbolic space. The algorithm is capable of handling surfaces of any genus and can work with singular triangulations. It also introduces a method to compute conjugate harmonic maps starting from a discrete harmonic map, which is used to generate conjugate minimal surfaces. This process preserves symmetry properties of the boundary curves and can yield unstable solutions for free boundary problems. The algorithm is based on minimizing the Dirichlet energy, which is a measure of the distortion of a map. It involves iteratively minimizing the energy by adjusting points in the image space and by varying points in the parameter domain. The algorithm is designed to handle both fixed and free boundary conditions, and it respects the symmetry properties of the boundary curves. The paper also discusses the conjugation algorithm, which transforms a discrete harmonic map into a conjugate harmonic map. This process is particularly useful for generating conjugate minimal surfaces from minimal surfaces. The conjugation algorithm preserves the symmetry properties of the original surface and is based on the geometric properties of the discrete surface. The algorithm is tested on various examples, including surfaces with different topologies and boundary conditions. The results demonstrate the effectiveness of the algorithm in computing discrete minimal surfaces and their conjugates, even in the presence of singularities and complex boundary conditions. The algorithm is implemented using the mathematical programming environment Grape, which is designed for problems in differential geometry and continuum mechanics.This paper presents a new algorithm for computing stable discrete minimal surfaces bounded by fixed or free boundary curves in 3D space, as well as in the 3-dimensional sphere and hyperbolic space. The algorithm is capable of handling surfaces of any genus and can work with singular triangulations. It also introduces a method to compute conjugate harmonic maps starting from a discrete harmonic map, which is used to generate conjugate minimal surfaces. This process preserves symmetry properties of the boundary curves and can yield unstable solutions for free boundary problems. The algorithm is based on minimizing the Dirichlet energy, which is a measure of the distortion of a map. It involves iteratively minimizing the energy by adjusting points in the image space and by varying points in the parameter domain. The algorithm is designed to handle both fixed and free boundary conditions, and it respects the symmetry properties of the boundary curves. The paper also discusses the conjugation algorithm, which transforms a discrete harmonic map into a conjugate harmonic map. This process is particularly useful for generating conjugate minimal surfaces from minimal surfaces. The conjugation algorithm preserves the symmetry properties of the original surface and is based on the geometric properties of the discrete surface. The algorithm is tested on various examples, including surfaces with different topologies and boundary conditions. The results demonstrate the effectiveness of the algorithm in computing discrete minimal surfaces and their conjugates, even in the presence of singularities and complex boundary conditions. The algorithm is implemented using the mathematical programming environment Grape, which is designed for problems in differential geometry and continuum mechanics.