u0 = 9.274E-24; % J/Tesla
a = 1/2*6.835E9; % Hz
h = 6.646E-34;
gI = 0.0009954;
gJ = 2.002331;

Iz = 1/2*[1 0; 0 -1]; Ip = [0 1; 0 0]; Im = [0 0; 1 0];
Jz = 1/2*[1 0; 0 -1]; Jp = [0 1; 0 0]; Jm = [0 0; 1 0];

IJ = 1/2*(kron(Ip,Jm)+kron(Im,Jp))+kron(Iz,Jz);

Bz = 0; % tesla
H = a*h*IJ + gJ*u0*kron(eye(2),Jz)*Bz - gI*u0*kron(Iz,eye(2))*Bz;