## 4 Numerical Example: Transient Fin Problem  
We implemented the method in 2D User Element Subroutine (UEL) within the commercial finite element analysis software Abaqus to demonstrate the performance of the Transient Thermal ZFEM method. Abaqus cannot solve systems with complex variables directly because its built-in equation solver is made for real-valued components. As a result, the complex-variable right-hand side of the element (\(RHS^{\epsilon}\)), and the complex-variable coefficient matrix for the element (\(A^{\epsilon}\)) are transformed to the Cauchy-Riemann form of all real numbers before returning the results to Abaqus. Then, the global system of equations is assembled internally by Abaqus and solved for each time increment. This process is presented in the flowchart of Figure 3, and the full version of the UEL is presented in the Supplementary Material. The method was verified by analyzing the heat transfer behavior for the thin fin shown in Figure 1a. The geometry of the fin is characterized by the length \(b\), height \(\delta\), and thickness \(L\), and it is made from a material with conductivity \(k\), density \(\rho\), and specific heat \(C_{P}\). The results computed using ZFEM were compared with the analytical solution using symbolic differentiation and Abaqus simulation with 2D DC2D4 elements and finite differentiation. The errors were measured using the Normalized Root Mean Square Error (NRMSE) metric as described in Equation (18). In this equation, \(o\) corresponds to the order of the sensitivity of \(g\), \(g^{(o)}\) is the reference results, and \(\#inc\) is the total number of sampling points which is coincident with the total number of increments in the transient thermal problem. \[\text{NRMSE}=\frac{1}{\text{max}\Big{(}g^{(o)}\Big{)}-\text{min}\Big{(}g^{(o)} \Big{)}}\sqrt{\frac{\sum\Big{(}g^{(o)}-g^{(o)}\Big{)}^{2}}{\#inc}} \tag{18}\]  
The analytical solution for the temperature as a function of time for the fin problem was obtained by considering the dimensionless parameters for the position \(X\) in Equation (19), the time \(\tau\) in Equation (20), \(N^{2}\) in Equation (21), \(\lambda_{n}\) in Equation (22), the temperature \(\theta\) in Equation (23), and by adding the steady-state \(\theta_{ss}(X)\) and the transient solutions \(\theta_{\tau}(X,\tau)\) as shown in Equations (24) and (25), respectively [28]. For this solution, it was also assumed  
Figure 3: Flowchart of Transient Thermal ZFEM and Abaqus User Element Subroutine. All complex operations are performed inside the UEL and the \(A\) matrix and the \(RHS\) vector are decomposed into Cauchy-Riemann form before exiting the UEL. that the axis \(x\) and \(X\) were aligned with the length of the fin, and the fin had an adiabatic tip. \[X=\frac{x}{b} \tag{19}\]  
\[\tau=\frac{tk}{b^{2}\rho C_{P}} \tag{20}\]  
\[N^{2}=\frac{2h_{U}b^{2}}{k\delta} \tag{21}\]  
\[\lambda_{n}=\frac{2n-1}{2}\tau \tag{22}\]  
\[\theta=\frac{T-T_{\infty}}{T_{W}-T_{\infty}}=\theta_{ss}(X)+\theta_{\tau}(X,\tau) \tag{23}\]  
\[\theta_{ss}(X)=\frac{\cosh(NX)}{\cosh(N)} \tag{24}\]  
\[\theta_{\tau}(X,\tau)=2\sum_{n=1}^{\infty}(-1)^{n+1}\bigg{(}\frac{\theta_{0}}{ \lambda_{n}}-\frac{\lambda_{n}^{2}}{N^{2}+\lambda_{n}^{2}}\bigg{)}cos(\lambda _{n}X)e^{-(N^{2}+\lambda_{n}^{2})\tau} \tag{25}\]  
In this problem, the fin was in thermal equilibrium with the environment at \(t=0\) s, meaning that the room temperature (\(T_{\infty}\)) was equal to the initial temperature (\(T_{0}\)). Then, a prescribed temperature \(T_{w}\) was applied at the base of the fin, i.e., \(x=b\). As soon as the fin started to heat, energy dissipation occurred due to the film conditions in the top, bottom, back, and front surfaces determined by \(T_{\infty}\) and the global heat transfer coefficient \(h_{U}\). In the case of ZFEM and Abaqus built-in models, combined convection and radiation boundary conditions were applied on the top and bottom edges, and \(T_{w}\) was applied to the left edge as depicted in Figure 1a.