We need to find eigenvalues and eigenstates of the hamiltonian H = A I.S + g μ B.S. At very  low field the eigenstates of this Hamiltonian are states defined by definite total angular momentum F = I+S, |F,(Subscript[m, F])>. At higher fields there is mixing between different states, but one can still describe them by F and (Subscript[m, F]) quantum numbers from which they evolved from low field.

For I =  3/2  and S = 
          1/2  there are a total of 8 states | F,(Subscript[m, F])>. Here is list of Kat states

i = 3/2;
s = 1/2;

Klist = {K[2, -2], K[2, -1], K[1, -1], K[2, 0], K[1, 0], 
        K[2, 1], K[1, 1], K[2, 2]}

We need to evaluate how the Hamiltonian acts on the states. We will use Non-commutative multiplication ** to indicate an operator acting on a state.

IS ** 
        K[F_, m_] ^= (F (F + 1) - i (i + 1) - s (s + 1))/
          2 K[F, m]

Here we use a slightly tricky definition of an upvalue (similar to object oriented programming), that associates this definition with the operator IS.

Sz ** 
        K[2, m_] ^= (m K[2, m] + Sqrt[2^2 - m^2] K[1, m])/(2  i + 
            1)
Sz ** 
        K[1, m_] ^= (Sqrt[2^2 - m^2] K[2, m] - m K[1, m])/(2  i + 
            1)

This operator requires some Clebsch-Gordon coefficients

HK = A Thread[IS ** Klist] +  gm Bz Thread[Sz ** Klist]

We can pick the z axis to point along B field

Tomat[e_] := 
      Table[Coefficient[e, Klist[[j]]], {j, 1, 8}]

This function will help build a matix of the Hamiltonian by acting with a "Bra" on the left hand side

HM = Tomat[HK]

MatrixForm[HM]

The form is block-diagonal with a series of 2x2 matrixes, so this can be handled analytically in principle.

EV = Eigenvalues[HM]

Table[
        Solve[D[EV[[j]] - EV[[k]], Bz] == 0, Bz], {j, 1, 
          8}, {k, 1, j}] // MatrixForm

General way of finding field-insensitive transitions

Table[
          Solve[D[EV[[j]] - EV[[k]], Bz] == 0, Bz], {j, 1, 
            8}, {k, 1, j}] /. {A -> 10^9, gm -> 2.8 10^6} // 
      MatrixForm

Plot[
      Evaluate[EV /. {A -> 10^9, gm -> 2.8 10^6}], {Bz, 0, 
        2000}]

We get a usual Breit-Rabi diagram

ES = Eigensystem[HM]

This makes it easy to see  which eigenvalue corresponds to which eigenstate. We care about allowed transitions with m changing by 0 or 1.

Plot[
      Evaluate[{EV[[1]] - EV[[5]], EV[[1]] - EV[[6]], 
            EV[[2]] - EV[[7]], 
            EV[[2]] - EV[[8]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

No interesting transitions among the ones originating in |2,2> or |2,-2> states

Plot[
      Evaluate[{EV[[3]] - EV[[4]], EV[[3]] - EV[[5]], 
            EV[[3]] - EV[[6]], 
            EV[[3]] - EV[[7]], 
            EV[[3]] - EV[[8]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

Here we got one among transitions out of the |1,0> state.

Plot[
      Evaluate[{EV[[3]] - EV[[6]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

Solve[
      D[Evaluate[{EV[[3]] - EV[[6]]} /. {A -> 
                  10^9, gm -> 2.8 10^6}], Bz] == 0, 
      Bz]

ES[[2, 3]] /. {A -> 10^9, gm -> 2.8 10^6, 
        Bz -> 191.392}

Mathamatica given  eigenvectors which are not normalized

Plot[
      Evaluate[{EV[[3]] - EV[[5]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 4000}]

Another field insensitive transition at much higher field

Plot[
      Evaluate[{EV[[4]] - EV[[5]], EV[[4]] - EV[[6]], 
            EV[[4]] - EV[[7]], 
            EV[[4]] - EV[[8]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

Here we got one among transitions out of the |2,0>state.

Plot[
      Evaluate[{EV[[4]] - EV[[5]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

Solve[
      D[Evaluate[{EV[[4]] - EV[[5]]} /. {A -> 10^9, 
                gm -> 2.8 10^6}], Bz] == 0, Bz]

ES[[2, 4]] /. {A -> 10^9, gm -> 2.8 10^6, 
        Bz -> 191.392}

Plot[
      Evaluate[{EV[[5]] - EV[[6]], 
            EV[[7]] - EV[[8]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

Plot[
      Evaluate[{EV[[5]] - EV[[6]]} /. {A -> 10^9, 
            gm -> 2.8 10^6}], {Bz, 0, 1000}]

Here is one between |2,-1> and |1,-1> states

Solve[
      D[Evaluate[{EV[[5]] - EV[[6]]} /. {A -> 10^9, 
                gm -> 2.8 10^6}], Bz] == 0, Bz]

ES[[2, 5]] /. {A -> 10^9, gm -> 2.8 10^6, 
        Bz -> 357.142}

So there is a total of 4  field-independent transitions at finite field, two at a field of 191 G and one at a field of 357 G and one at field of 2666 G