The problem of numerically integrating general elliptic differential equations in irregular two and three dimensional regions is discussed. The method used numerically computes a transformation of the given region into a rectangular region. The numerical coordinate transformation is determined by requiring that the components of the transformation satisfy inhomogeneous Laplace or more general equations. The transformation is then used to transform the different equation and the boundary conditions to the rectangular region. The boundary value problem in the rectangular region is integrated using one of the standard methods for general elliptic equations. The use of the existing software reduces the problem to analytically transforming the given differential equation and the Laplacian to general coordinate frame and then writing subroutines that will tabulate the coefficients of these differential equations using the tabulated coordinate transformation. This method has been successfully used in two dimensions so we are concerned with the three-dimensional extension of the existing codes where the major problem encountered is the volume of algebra and coding required to complete the method. To overcome these difficulties, a symbol manipulation program in VAXIMA is written that has an input the formula for the given differential equation in some natural coordinates and has as output the required FORTRAN subroutines. Because of the complexity of the resulting code, code validation was performed by systematic truncation error testing. The paper concludes with a discussion of the problems encountered in using a symbol manipulator to write large FORTRAN codes.