This simulation study evaluated the effect of three different inlay materials (composite, glass ceramic, zirconia), cavity design parameters (isthmus width and depth) and different force vectors and magnitudes on the stress distribution within mesio-occlusal-distal (MOD) inlays and the remaining enamel and dentin. The mechanical performance of inlays was evaluated using 3-D finite element analysis (FEA) method. Three different restoration materials and hard tissues of the restored tooth with different cavity depth (2–5 mm) and width (2–4 mm) were exposed to occlusal loading with different magnitudes from 10 to 130 kg at varying angles between 0° and 15°. The maximum von Mises stresses were calculated for the inlays, tooth structure and bonded surfaces. Response Surface Optimization method was implemented into the finite element software package in order to design cavity shapes with more favourable interfacial stresses for bonded restorations under occlusal loading. Teeth restored with resin composite exhibited the highest von Mises Stress, followed by glass ceramic and zirconia. The increase in isthmus width decreased interfacial shear stresses in zirconia MOD inlay but the increase in cavity depth did not change the stress levels for all three materials. According to mechanical safety factor, inlay and tooth structure remained within the mechanical limits in three parameters (material, magnitude of force, cavity shape) but negatively affected by the force vector.