We have developed a flexible forward modeling program for fault-bend folding based on a velocity description of deformation. The program can handle multiple fault bends of any geometry (e. g. fault bends not stepping up from a detachment) provided they are not wedge geometries, and variable active axial surface orientation with corresponding varying slip ratios. For fault shapes in which the fundamental fault-bend fold equations are valid the program can produce hanging wall fault-bend folds that conserve cross-sectional area, layer thickness, and line length. For fault shapes in which there is no solution (e. g. a ramp > 30 deg. flattening to a detachment) the program can solve for an axial surface orientation that minimizes changes in layer thickness and line length while still conserving cross-sectional area, thus providing a reasonable end-member solution that minimizes non-flexural slip strain. In addition, because the active axial surface orientation can be independently defined from fault shape, the program can be used in an inclined shear fashion (as is often desirable in modeling extensional fault-bend folding). The program stores the entire deformation history of the hanging wall for all time steps for all particle locations in a hierarchical array, allowing for the exploration of any parameter to be analyzed through time (such as change in line length of a layer, or variation in slip of a particle). This program allows for the easy exploration of fault-bend fold geometries that can be compared to suspected real fault-bend folds, and perhaps more importantly provides a framework for future development of an inverse modeling program for fault-bend folding.