def METHOD_NAME(kbu, f_obs):
  p = flex.double([10, 20])
  kbu.update(b_sols = p)
  tg = bulk_solvent.ls_kbp_sol_u_star(
    f_model     = kbu.data,
    f_obs       = f_obs.data(),
    scale       = 1.0,
    kb_sol_grad = True,
    p_sol_grad  = False,
    u_star_grad = True,
    kb_sol_curv = True,
    p_sol_curv  = False)
  g_b_sols_anal = list(tg.grad_b_sols())
  c_b_sols_anal = list(tg.curv_b_sols())
  print("C anal b:", c_b_sols_anal)
  #
  g_b_sols_fd = []
  c_b_sols_fd = []
  t0 = tg.target()
  for shift in [flex.double([1.e-3,0]),
                flex.double([0,1.e-3])]:
    # plus shift
    kbu.update(b_sols = p+shift)
    tg = bulk_solvent.ls_kbp_sol_u_star(
      f_model     = kbu.data,
      f_obs       = f_obs.data(),
      scale       = 1.0,
      kb_sol_grad = True,
      p_sol_grad  = False,
      u_star_grad = True,
      kb_sol_curv = True,
      p_sol_curv  = False)
    t1 = tg.target()
    # minus shift
    kbu.update(b_sols = p-shift)
    tg = bulk_solvent.ls_kbp_sol_u_star(
      f_model     = kbu.data,
      f_obs       = f_obs.data(),
      scale       = 1.0,
      kb_sol_grad = True,
      p_sol_grad  = False,
      u_star_grad = True,
      kb_sol_curv = True,
      p_sol_curv  = False)
    t2 = tg.target()
    g_b_sols_fd.append( (t1-t2)/(2*1.e-3) )
    c_b_sols_fd.append( (t1+t2-2*t0)/1.e-3**2 *2)
    #
  print("C df:",c_b_sols_fd)
  assert approx_equal(g_b_sols_anal, g_b_sols_fd)