We present VAH, a (3+1)–dimensional simulation that evolves the far-from-equilibrium quark-gluon plasma produced in ultrarelativistic heavy-ion collisions with anisotropic fluid dynamics. We solve the hydrodynamic equations on an Eulerian grid using the Kurganov–Tadmor algorithm in combination with a new adaptive Runge–Kutta method. Our numerical scheme allows us to start the simulation soon after the nuclear collision, largely avoiding the need to integrate it with a separate pre-equilibrium dynamics module. We test the code's performance by simulating on the Eulerian grid conformal and non-conformal Bjorken flow as well as conformal Gubser flow, whose (0+1)–dimensional solutions are precisely known. Finally, we compare non-conformal anisotropic hydrodynamics to second-order viscous hydrodynamics in central Pb+Pb collisions and find that the former's longitudinal flow profile responds more consistently to the fluid's gradients along the spacetime rapidity direction. Program summaryManuscript Title: Anisotropic fluid dynamical simulations of heavy-ion collisionsAuthors: Mike McNelis, Dennis Bazow, Ulrich HeinzProgram Title:VAHLicensing provisions: GPLv3Programming Language: C++Computer: Laptop, desktop, clusterOperating System: GNU/Linux distributions, Mac OS XGlobal memory usage: 1.2 GB (for a 129×129×63 grid)Keywords: Ultrarelativistic heavy-ion collisions, quark-gluon plasma, relativistic hydrodynamics, computational fluid dynamicsClassification: 12 Gases and Fluids, 17 Nuclear physicsExternal routines/libraries: GNU Scientific Library (GSL)Nature of problem: Modeling the far-from-equilibrium dynamics of quark-gluon plasma produced in ultrarelativistic heavy-ion collisions.Solution method: Kurganov–Tadmor algorithm, adaptive stepsize methodRunning time: A (3+1)–d non-conformal anisotropic fluid dynamical simulation of a central Pb+Pb collision on a 129×129×63 grid takes about 530s for an Intel Xeon E5-2680 v4 multi-core processor with OpenMP acceleration (see Sec. 5.3).