<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def get_delta_plus_epsilon(nat, i_now, dr_now, dr_old):
 """ Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix """

    XIJ = get_XIJ(nat, i_now, dr_now, dr_old)
    YIJ = get_YIJ(nat, i_now, dr_old)

    YIJ_invert = array_inverse(YIJ)

    # Perform sum_k X_ik Y_jk^-1
    epsilon = np.sum(XIJ.reshape(-1,3,1,3)*YIJ_invert.reshape(-1,1,3,3), axis=3)

    return epsilon