How to solve a system of differential equations using ode45 and vary the magnetic field and plot the solutions as a function of this field?

1 view (last 30 days)
I'm trying to solve this system of differential equations, but I can only solve it if I keep B constant. I wanted to vary B then plot the solutions as a function of B. If anyone can help me I would be very grateful.
B=[0:0.1:10];
d0=0.2
wb=1
wd=wb-d0
wl=wb
db=0.18
dd=0.05
a=0.02
Ua=0.15
Ub=0.02
g1= 0.001
g2=g1
g3= 0.0001
g4=g3
%syms B theta mub ghx ghz gex gez d0 wb wd wl db dd a Ua Ub g1 g2 g3 g4
theta=pi/3
mub=0.0579; %# Bohr magneton
ghx=-0.35 ;
ghz=-2.2;
gex=-0.65;
gez=-0.8;
bp=mub*B.*sin(theta)*(gez+ghz)*0.5 ; %bp
bm=mub*B.*sin(theta)*(gez-ghz)*0.5; %bm
be=mub*B.*cos(theta)*gex*0.5; %be
bh=mub*B.*cos(theta)*ghx*0.5; %bh
dRdt=@(t,R)[R(7)*g1 + R(13)*g2 + R(19)*g3 + R(25)*g4 + R(3)*Ua*1i - R(11)*Ua*1i + R(4)*Ub*1i - R(16)*Ub*1i
- R(12)*Ua*1i - R(17)*Ub*1i + (R(3)*db*1i)/2 - (R(2)*g1)/2 + R(2)*(a*B.^2 + (mub*sin(theta)*(gez + ghz)*B)/2 + wb - wl)*1i + (B.*R(4)*gex*mub*cos(theta)*1i)/2 + (B.*R(5)*ghx*mub*cos(theta)*1i)/2
R(1)*Ua*1i - R(13)*Ua*1i - R(18)*Ub*1i + (R(2)*db*1i)/2 - (R(3)*g2)/2 + R(3)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(5)*gex*mub*cos(theta)*1i)/2 + (B.*R(4)*ghx*mub*cos(theta)*1i)/2
- R(14)*Ua*1i + R(1)*Ub*1i - R(19)*Ub*1i + (R(5)*dd*1i)/2 - (R(4)*g3)/2 + R(4)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + (B.*R(2)*gex*mub*cos(theta)*1i)/2 + (B.*R(3)*ghx*mub*cos(theta)*1i)/2
- R(15)*Ua*1i - R(20)*Ub*1i + (R(4)*dd*1i)/2 - (R(5)*g4)/2 + R(5)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + (B.*R(3)*gex*mub*cos(theta)*1i)/2 + (B.*R(2)*ghx*mub*cos(theta)*1i)/2
R(8)*Ua*1i + R(9)*Ub*1i - (R(11)*db*1i)/2 - (R(6)*g1)/2 - R(6)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(16)*gex*mub*cos(theta)*1i)/2 - (B.*R(21)*ghx*mub*cos(theta)*1i)/2
(R(8)*db*1i)/2 - (R(12)*db*1i)/2 - R(7)*g1 + (B.*R(9)*gex*mub*cos(theta)*1i)/2 - (B.*R(17)*gex*mub*cos(theta)*1i)/2 + (B.*R(10)*ghx*mub*cos(theta)*1i)/2 - (B.*R(22)*ghx*mub*cos(theta)*1i)/2
R(6)*Ua*1i + (R(7)*db*1i)/2 - (R(13)*db*1i)/2 - (R(8)*g1)/2 - (R(8)*g2)/2 + R(8)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - R(8)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(10)*gex*mub*cos(theta)*1i)/2 - (B*R(18)*gex*mub*cos(theta)*1i)/2 + (B.*R(9)*ghx*mub*cos(theta)*1i)/2 - (B.*R(23)*ghx*mub*cos(theta)*1i)/2
R(6)*Ub*1i - (R(14)*db*1i)/2 + (R(10)*dd*1i)/2 - (R(9)*g1)/2 - (R(9)*g3)/2 + R(9)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(9)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(7)*gex*mub*cos(theta)*1i)/2 - (B.*R(19)*gex*mub*cos(theta)*1i)/2 + (B.*R(8)*ghx*mub*cos(theta)*1i)/2 - (B.*R(24)*ghx*mub*cos(theta)*1i)/2
- (R(15)*db*1i)/2 + (R(9)*dd*1i)/2 - (R(10)*g1)/2 - (R(10)*g4)/2 + R(10)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(10)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(8)*gex*mub*cos(theta)*1i)/2 - (B.*R(20)*gex*mub*cos(theta)*1i)/2 + (B.*R(7)*ghx*mub*cos(theta)*1i)/2 - (B.*R(25)*ghx*mub*cos(theta)*1i)/2
- R(1)*Ua*1i + R(13)*Ua*1i + R(14)*Ub*1i - (R(6)*db*1i)/2 - (R(11)*g2)/2 - R(11)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(21)*gex*mub*cos(theta)*1i)/2 - (B.*R(16)*ghx*mub*cos(theta)*1i)/2
- R(2)*Ua*1i - (R(7)*db*1i)/2 + (R(13)*db*1i)/2 - (R(12)*g1)/2 - (R(12)*g2)/2 - R(12)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + R(12)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(14)*gex*mub*cos(theta)*1i)/2 - (B.*R(22)*gex*mub*cos(theta)*1i)/2 + (B.*R(15)*ghx*mub*cos(theta)*1i)/2 - (B.*R(17)*ghx*mub*cos(theta)*1i)/2
- R(3)*Ua*1i + R(11)*Ua*1i - (R(8)*db*1i)/2 + (R(12)*db*1i)/2 - R(13)*g2 + (B.*R(15)*gex*mub*cos(theta)*1i)/2 - (B.*R(23)*gex*mub*cos(theta)*1i)/2 + (B.*R(14)*ghx*mub*cos(theta)*1i)/2 - (B.*R(18)*ghx*mub*cos(theta)*1i)/2
- R(4)*Ua*1i + R(11)*Ub*1i - (R(9)*db*1i)/2 + (R(15)*dd*1i)/2 - (R(14)*g2)/2 - (R(14)*g3)/2 + R(14)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(14)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(12)*gex*mub*cos(theta)*1i)/2 - (B.*R(24)*gex*mub*cos(theta)*1i)/2 + (B.*R(13)*ghx*mub*cos(theta)*1i)/2 - (B.*R(19)*ghx*mub*cos(theta)*1i)/2
- R(5)*Ua*1i - (R(10)*db*1i)/2 + (R(14)*dd*1i)/2 - (R(15)*g2)/2 - (R(15)*g4)/2 + R(15)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(15)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(13)*gex*mub*cos(theta)*1i)/2 - (B.*R(25)*gex*mub*cos(theta)*1i)/2 + (B.*R(12)*ghx*mub*cos(theta)*1i)/2 - (B.*R(20)*ghx*mub*cos(theta)*1i)/2
R(18)*Ua*1i - R(1)*Ub*1i + R(19)*Ub*1i - (R(21)*dd*1i)/2 - (R(16)*g3)/2 - R(16)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(6)*gex*mub*cos(theta)*1i)/2 - (B.*R(11)*ghx*mub*cos(theta)*1i)/2
- R(2)*Ub*1i + (R(18)*db*1i)/2 - (R(22)*dd*1i)/2 - (R(17)*g1)/2 - (R(17)*g3)/2 - R(17)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(17)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(7)*gex*mub*cos(theta)*1i)/2 + (B.*R(19)*gex*mub*cos(theta)*1i)/2 - (B.*R(12)*ghx*mub*cos(theta)*1i)/2 + (B.*R(20)*ghx*mub*cos(theta)*1i)/2
R(16)*Ua*1i - R(3)*Ub*1i + (R(17)*db*1i)/2 - (R(23)*dd*1i)/2 - (R(18)*g2)/2 - (R(18)*g3)/2 - R(18)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(18)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(8)*gex*mub*cos(theta)*1i)/2 + (B.*R(20)*gex*mub*cos(theta)*1i)/2 - (B.*R(13)*ghx*mub*cos(theta)*1i)/2 + (B.*R(19)*ghx*mub*cos(theta)*1i)/2
- R(4)*Ub*1i + R(16)*Ub*1i + (R(20)*dd*1i)/2 - (R(24)*dd*1i)/2 - R(19)*g3 - (B.*R(9)*gex*mub*cos(theta)*1i)/2 + (B.*R(17)*gex*mub*cos(theta)*1i)/2 - (B.*R(14)*ghx*mub*cos(theta)*1i)/2 + (B.*R(18)*ghx*mub*cos(theta)*1i)/2
- R(5)*Ub*1i + (R(19)*dd*1i)/2 - (R(25)*dd*1i)/2 - (R(20)*g3)/2 - (R(20)*g4)/2 - R(20)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(20)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(10)*gex*mub*cos(theta)*1i)/2 + (B.*R(18)*gex*mub*cos(theta)*1i)/2 - (B.*R(15)*ghx*mub*cos(theta)*1i)/2 + (B.*R(17)*ghx*mub*cos(theta)*1i)/2
R(23)*Ua*1i + R(24)*Ub*1i - (R(16)*dd*1i)/2 - (R(21)*g4)/2 - R(21)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(11)*gex*mub*cos(theta)*1i)/2 - (B.*R(6)*ghx*mub*cos(theta)*1i)/2
(R(23)*db*1i)/2 - (R(17)*dd*1i)/2 - (R(22)*g1)/2 - (R(22)*g4)/2 - R(22)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(22)*(a.*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(12)*gex*mub*cos(theta)*1i)/2 + (B.*R(24)*gex*mub*cos(theta)*1i)/2 - (B.*R(7)*ghx*mub*cos(theta)*1i)/2 + (B.*R(25)*ghx*mub*cos(theta)*1i)/2
R(21)*Ua*1i + (R(22)*db*1i)/2 - (R(18)*dd*1i)/2 - (R(23)*g2)/2 - (R(23)*g4)/2 - R(23)*(a.*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(23)*(a.*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(13)*gex*mub*cos(theta)*1i)/2 + (B.*R(25)*gex*mub*cos(theta)*1i)/2 - (B.*R(8)*ghx*mub*cos(theta)*1i)/2 + (B.*R(24)*ghx*mub*cos(theta)*1i)/2
R(21)*Ub*1i - (R(19)*dd*1i)/2 + (R(25)*dd*1i)/2 - (R(24)*g3)/2 - (R(24)*g4)/2 + R(24)*(a.*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(24)*(a.*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(14)*gex*mub*cos(theta)*1i)/2 + (B.*R(22)*gex*mub*cos(theta)*1i)/2 - (B.*R(9)*ghx*mub*cos(theta)*1i)/2 + (B.*R(23)*ghx*mub*cos(theta)*1i)/2
- (R(20)*dd*1i)/2 + (R(24)*dd*1i)/2 - R(25)*g4 - (B.*R(15)*gex*mub*cos(theta)*1i)/2 + (B.*R(23)*gex*mub*cos(theta)*1i)/2 - (B.*R(10)*ghx*mub*cos(theta)*1i)/2 + (B.*R(22)*ghx*mub*cos(theta)*1i)/2]
[t,R]=ode45(dRdt,[0 100],[1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0])
plot(B,R(:,1))

Accepted Answer

Star Strider
Star Strider on 15 May 2022
First, define ‘dRdt’ as:
dRdt=@(t,R,B)
then loop to substitute differeent values of ‘B’ and store the results.
B=[0:0.1:10];
B = linspace(0, 10, 200);
d0=0.2
d0 = 0.2000
wb=1
wb = 1
wd=wb-d0
wd = 0.8000
wl=wb
wl = 1
db=0.18
db = 0.1800
dd=0.05
dd = 0.0500
a=0.02
a = 0.0200
Ua=0.15
Ua = 0.1500
Ub=0.02
Ub = 0.0200
g1= 0.001
g1 = 1.0000e-03
g2=g1
g2 = 1.0000e-03
g3= 0.0001
g3 = 1.0000e-04
g4=g3
g4 = 1.0000e-04
%syms B theta mub ghx ghz gex gez d0 wb wd wl db dd a Ua Ub g1 g2 g3 g4
theta=pi/3
theta = 1.0472
mub=0.0579; %# Bohr magneton
ghx=-0.35 ;
ghz=-2.2;
gex=-0.65;
gez=-0.8;
bp=mub*B.*sin(theta)*(gez+ghz)*0.5 ; %bp
bm=mub*B.*sin(theta)*(gez-ghz)*0.5; %bm
be=mub*B.*cos(theta)*gex*0.5; %be
bh=mub*B.*cos(theta)*ghx*0.5; %bh
dRdt=@(t,R,B)[R(7)*g1 + R(13)*g2 + R(19)*g3 + R(25)*g4 + R(3)*Ua*1i - R(11)*Ua*1i + R(4)*Ub*1i - R(16)*Ub*1i
- R(12)*Ua*1i - R(17)*Ub*1i + (R(3)*db*1i)/2 - (R(2)*g1)/2 + R(2)*(a*B.^2 + (mub*sin(theta)*(gez + ghz)*B)/2 + wb - wl)*1i + (B.*R(4)*gex*mub*cos(theta)*1i)/2 + (B.*R(5)*ghx*mub*cos(theta)*1i)/2
R(1)*Ua*1i - R(13)*Ua*1i - R(18)*Ub*1i + (R(2)*db*1i)/2 - (R(3)*g2)/2 + R(3)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(5)*gex*mub*cos(theta)*1i)/2 + (B.*R(4)*ghx*mub*cos(theta)*1i)/2
- R(14)*Ua*1i + R(1)*Ub*1i - R(19)*Ub*1i + (R(5)*dd*1i)/2 - (R(4)*g3)/2 + R(4)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + (B.*R(2)*gex*mub*cos(theta)*1i)/2 + (B.*R(3)*ghx*mub*cos(theta)*1i)/2
- R(15)*Ua*1i - R(20)*Ub*1i + (R(4)*dd*1i)/2 - (R(5)*g4)/2 + R(5)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + (B.*R(3)*gex*mub*cos(theta)*1i)/2 + (B.*R(2)*ghx*mub*cos(theta)*1i)/2
R(8)*Ua*1i + R(9)*Ub*1i - (R(11)*db*1i)/2 - (R(6)*g1)/2 - R(6)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(16)*gex*mub*cos(theta)*1i)/2 - (B.*R(21)*ghx*mub*cos(theta)*1i)/2
(R(8)*db*1i)/2 - (R(12)*db*1i)/2 - R(7)*g1 + (B.*R(9)*gex*mub*cos(theta)*1i)/2 - (B.*R(17)*gex*mub*cos(theta)*1i)/2 + (B.*R(10)*ghx*mub*cos(theta)*1i)/2 - (B.*R(22)*ghx*mub*cos(theta)*1i)/2
R(6)*Ua*1i + (R(7)*db*1i)/2 - (R(13)*db*1i)/2 - (R(8)*g1)/2 - (R(8)*g2)/2 + R(8)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - R(8)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(10)*gex*mub*cos(theta)*1i)/2 - (B*R(18)*gex*mub*cos(theta)*1i)/2 + (B.*R(9)*ghx*mub*cos(theta)*1i)/2 - (B.*R(23)*ghx*mub*cos(theta)*1i)/2
R(6)*Ub*1i - (R(14)*db*1i)/2 + (R(10)*dd*1i)/2 - (R(9)*g1)/2 - (R(9)*g3)/2 + R(9)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(9)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(7)*gex*mub*cos(theta)*1i)/2 - (B.*R(19)*gex*mub*cos(theta)*1i)/2 + (B.*R(8)*ghx*mub*cos(theta)*1i)/2 - (B.*R(24)*ghx*mub*cos(theta)*1i)/2
- (R(15)*db*1i)/2 + (R(9)*dd*1i)/2 - (R(10)*g1)/2 - (R(10)*g4)/2 + R(10)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(10)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(8)*gex*mub*cos(theta)*1i)/2 - (B.*R(20)*gex*mub*cos(theta)*1i)/2 + (B.*R(7)*ghx*mub*cos(theta)*1i)/2 - (B.*R(25)*ghx*mub*cos(theta)*1i)/2
- R(1)*Ua*1i + R(13)*Ua*1i + R(14)*Ub*1i - (R(6)*db*1i)/2 - (R(11)*g2)/2 - R(11)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(21)*gex*mub*cos(theta)*1i)/2 - (B.*R(16)*ghx*mub*cos(theta)*1i)/2
- R(2)*Ua*1i - (R(7)*db*1i)/2 + (R(13)*db*1i)/2 - (R(12)*g1)/2 - (R(12)*g2)/2 - R(12)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + R(12)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(14)*gex*mub*cos(theta)*1i)/2 - (B.*R(22)*gex*mub*cos(theta)*1i)/2 + (B.*R(15)*ghx*mub*cos(theta)*1i)/2 - (B.*R(17)*ghx*mub*cos(theta)*1i)/2
- R(3)*Ua*1i + R(11)*Ua*1i - (R(8)*db*1i)/2 + (R(12)*db*1i)/2 - R(13)*g2 + (B.*R(15)*gex*mub*cos(theta)*1i)/2 - (B.*R(23)*gex*mub*cos(theta)*1i)/2 + (B.*R(14)*ghx*mub*cos(theta)*1i)/2 - (B.*R(18)*ghx*mub*cos(theta)*1i)/2
- R(4)*Ua*1i + R(11)*Ub*1i - (R(9)*db*1i)/2 + (R(15)*dd*1i)/2 - (R(14)*g2)/2 - (R(14)*g3)/2 + R(14)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(14)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(12)*gex*mub*cos(theta)*1i)/2 - (B.*R(24)*gex*mub*cos(theta)*1i)/2 + (B.*R(13)*ghx*mub*cos(theta)*1i)/2 - (B.*R(19)*ghx*mub*cos(theta)*1i)/2
- R(5)*Ua*1i - (R(10)*db*1i)/2 + (R(14)*dd*1i)/2 - (R(15)*g2)/2 - (R(15)*g4)/2 + R(15)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(15)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(13)*gex*mub*cos(theta)*1i)/2 - (B.*R(25)*gex*mub*cos(theta)*1i)/2 + (B.*R(12)*ghx*mub*cos(theta)*1i)/2 - (B.*R(20)*ghx*mub*cos(theta)*1i)/2
R(18)*Ua*1i - R(1)*Ub*1i + R(19)*Ub*1i - (R(21)*dd*1i)/2 - (R(16)*g3)/2 - R(16)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(6)*gex*mub*cos(theta)*1i)/2 - (B.*R(11)*ghx*mub*cos(theta)*1i)/2
- R(2)*Ub*1i + (R(18)*db*1i)/2 - (R(22)*dd*1i)/2 - (R(17)*g1)/2 - (R(17)*g3)/2 - R(17)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(17)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(7)*gex*mub*cos(theta)*1i)/2 + (B.*R(19)*gex*mub*cos(theta)*1i)/2 - (B.*R(12)*ghx*mub*cos(theta)*1i)/2 + (B.*R(20)*ghx*mub*cos(theta)*1i)/2
R(16)*Ua*1i - R(3)*Ub*1i + (R(17)*db*1i)/2 - (R(23)*dd*1i)/2 - (R(18)*g2)/2 - (R(18)*g3)/2 - R(18)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(18)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(8)*gex*mub*cos(theta)*1i)/2 + (B.*R(20)*gex*mub*cos(theta)*1i)/2 - (B.*R(13)*ghx*mub*cos(theta)*1i)/2 + (B.*R(19)*ghx*mub*cos(theta)*1i)/2
- R(4)*Ub*1i + R(16)*Ub*1i + (R(20)*dd*1i)/2 - (R(24)*dd*1i)/2 - R(19)*g3 - (B.*R(9)*gex*mub*cos(theta)*1i)/2 + (B.*R(17)*gex*mub*cos(theta)*1i)/2 - (B.*R(14)*ghx*mub*cos(theta)*1i)/2 + (B.*R(18)*ghx*mub*cos(theta)*1i)/2
- R(5)*Ub*1i + (R(19)*dd*1i)/2 - (R(25)*dd*1i)/2 - (R(20)*g3)/2 - (R(20)*g4)/2 - R(20)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(20)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(10)*gex*mub*cos(theta)*1i)/2 + (B.*R(18)*gex*mub*cos(theta)*1i)/2 - (B.*R(15)*ghx*mub*cos(theta)*1i)/2 + (B.*R(17)*ghx*mub*cos(theta)*1i)/2
R(23)*Ua*1i + R(24)*Ub*1i - (R(16)*dd*1i)/2 - (R(21)*g4)/2 - R(21)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(11)*gex*mub*cos(theta)*1i)/2 - (B.*R(6)*ghx*mub*cos(theta)*1i)/2
(R(23)*db*1i)/2 - (R(17)*dd*1i)/2 - (R(22)*g1)/2 - (R(22)*g4)/2 - R(22)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(22)*(a.*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(12)*gex*mub*cos(theta)*1i)/2 + (B.*R(24)*gex*mub*cos(theta)*1i)/2 - (B.*R(7)*ghx*mub*cos(theta)*1i)/2 + (B.*R(25)*ghx*mub*cos(theta)*1i)/2
R(21)*Ua*1i + (R(22)*db*1i)/2 - (R(18)*dd*1i)/2 - (R(23)*g2)/2 - (R(23)*g4)/2 - R(23)*(a.*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(23)*(a.*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(13)*gex*mub*cos(theta)*1i)/2 + (B.*R(25)*gex*mub*cos(theta)*1i)/2 - (B.*R(8)*ghx*mub*cos(theta)*1i)/2 + (B.*R(24)*ghx*mub*cos(theta)*1i)/2
R(21)*Ub*1i - (R(19)*dd*1i)/2 + (R(25)*dd*1i)/2 - (R(24)*g3)/2 - (R(24)*g4)/2 + R(24)*(a.*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(24)*(a.*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(14)*gex*mub*cos(theta)*1i)/2 + (B.*R(22)*gex*mub*cos(theta)*1i)/2 - (B.*R(9)*ghx*mub*cos(theta)*1i)/2 + (B.*R(23)*ghx*mub*cos(theta)*1i)/2
- (R(20)*dd*1i)/2 + (R(24)*dd*1i)/2 - R(25)*g4 - (B.*R(15)*gex*mub*cos(theta)*1i)/2 + (B.*R(23)*gex*mub*cos(theta)*1i)/2 - (B.*R(10)*ghx*mub*cos(theta)*1i)/2 + (B.*R(22)*ghx*mub*cos(theta)*1i)/2]
dRdt = function_handle with value:
@(t,R,B)[R(7)*g1+R(13)*g2+R(19)*g3+R(25)*g4+R(3)*Ua*1i-R(11)*Ua*1i+R(4)*Ub*1i-R(16)*Ub*1i;-R(12)*Ua*1i-R(17)*Ub*1i+(R(3)*db*1i)/2-(R(2)*g1)/2+R(2)*(a*B.^2+(mub*sin(theta)*(gez+ghz)*B)/2+wb-wl)*1i+(B.*R(4)*gex*mub*cos(theta)*1i)/2+(B.*R(5)*ghx*mub*cos(theta)*1i)/2;R(1)*Ua*1i-R(13)*Ua*1i-R(18)*Ub*1i+(R(2)*db*1i)/2-(R(3)*g2)/2+R(3)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(5)*gex*mub*cos(theta)*1i)/2+(B.*R(4)*ghx*mub*cos(theta)*1i)/2;-R(14)*Ua*1i+R(1)*Ub*1i-R(19)*Ub*1i+(R(5)*dd*1i)/2-(R(4)*g3)/2+R(4)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+(B.*R(2)*gex*mub*cos(theta)*1i)/2+(B.*R(3)*ghx*mub*cos(theta)*1i)/2;-R(15)*Ua*1i-R(20)*Ub*1i+(R(4)*dd*1i)/2-(R(5)*g4)/2+R(5)*(a*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+(B.*R(3)*gex*mub*cos(theta)*1i)/2+(B.*R(2)*ghx*mub*cos(theta)*1i)/2;R(8)*Ua*1i+R(9)*Ub*1i-(R(11)*db*1i)/2-(R(6)*g1)/2-R(6)*(a*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-(B.*R(16)*gex*mub*cos(theta)*1i)/2-(B.*R(21)*ghx*mub*cos(theta)*1i)/2;(R(8)*db*1i)/2-(R(12)*db*1i)/2-R(7)*g1+(B.*R(9)*gex*mub*cos(theta)*1i)/2-(B.*R(17)*gex*mub*cos(theta)*1i)/2+(B.*R(10)*ghx*mub*cos(theta)*1i)/2-(B.*R(22)*ghx*mub*cos(theta)*1i)/2;R(6)*Ua*1i+(R(7)*db*1i)/2-(R(13)*db*1i)/2-(R(8)*g1)/2-(R(8)*g2)/2+R(8)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-R(8)*(a*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(10)*gex*mub*cos(theta)*1i)/2-(B*R(18)*gex*mub*cos(theta)*1i)/2+(B.*R(9)*ghx*mub*cos(theta)*1i)/2-(B.*R(23)*ghx*mub*cos(theta)*1i)/2;R(6)*Ub*1i-(R(14)*db*1i)/2+(R(10)*dd*1i)/2-(R(9)*g1)/2-(R(9)*g3)/2+R(9)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-R(9)*(a*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(7)*gex*mub*cos(theta)*1i)/2-(B.*R(19)*gex*mub*cos(theta)*1i)/2+(B.*R(8)*ghx*mub*cos(theta)*1i)/2-(B.*R(24)*ghx*mub*cos(theta)*1i)/2;-(R(15)*db*1i)/2+(R(9)*dd*1i)/2-(R(10)*g1)/2-(R(10)*g4)/2+R(10)*(a*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-R(10)*(a*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(8)*gex*mub*cos(theta)*1i)/2-(B.*R(20)*gex*mub*cos(theta)*1i)/2+(B.*R(7)*ghx*mub*cos(theta)*1i)/2-(B.*R(25)*ghx*mub*cos(theta)*1i)/2;-R(1)*Ua*1i+R(13)*Ua*1i+R(14)*Ub*1i-(R(6)*db*1i)/2-(R(11)*g2)/2-R(11)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-(B.*R(21)*gex*mub*cos(theta)*1i)/2-(B.*R(16)*ghx*mub*cos(theta)*1i)/2;-R(2)*Ua*1i-(R(7)*db*1i)/2+(R(13)*db*1i)/2-(R(12)*g1)/2-(R(12)*g2)/2-R(12)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+R(12)*(a*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(14)*gex*mub*cos(theta)*1i)/2-(B.*R(22)*gex*mub*cos(theta)*1i)/2+(B.*R(15)*ghx*mub*cos(theta)*1i)/2-(B.*R(17)*ghx*mub*cos(theta)*1i)/2;-R(3)*Ua*1i+R(11)*Ua*1i-(R(8)*db*1i)/2+(R(12)*db*1i)/2-R(13)*g2+(B.*R(15)*gex*mub*cos(theta)*1i)/2-(B.*R(23)*gex*mub*cos(theta)*1i)/2+(B.*R(14)*ghx*mub*cos(theta)*1i)/2-(B.*R(18)*ghx*mub*cos(theta)*1i)/2;-R(4)*Ua*1i+R(11)*Ub*1i-(R(9)*db*1i)/2+(R(15)*dd*1i)/2-(R(14)*g2)/2-(R(14)*g3)/2+R(14)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-R(14)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(12)*gex*mub*cos(theta)*1i)/2-(B.*R(24)*gex*mub*cos(theta)*1i)/2+(B.*R(13)*ghx*mub*cos(theta)*1i)/2-(B.*R(19)*ghx*mub*cos(theta)*1i)/2;-R(5)*Ua*1i-(R(10)*db*1i)/2+(R(14)*dd*1i)/2-(R(15)*g2)/2-(R(15)*g4)/2+R(15)*(a*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-R(15)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i+(B.*R(13)*gex*mub*cos(theta)*1i)/2-(B.*R(25)*gex*mub*cos(theta)*1i)/2+(B.*R(12)*ghx*mub*cos(theta)*1i)/2-(B.*R(20)*ghx*mub*cos(theta)*1i)/2;R(18)*Ua*1i-R(1)*Ub*1i+R(19)*Ub*1i-(R(21)*dd*1i)/2-(R(16)*g3)/2-R(16)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-(B.*R(6)*gex*mub*cos(theta)*1i)/2-(B.*R(11)*ghx*mub*cos(theta)*1i)/2;-R(2)*Ub*1i+(R(18)*db*1i)/2-(R(22)*dd*1i)/2-(R(17)*g1)/2-(R(17)*g3)/2-R(17)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+R(17)*(a*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-(B.*R(7)*gex*mub*cos(theta)*1i)/2+(B.*R(19)*gex*mub*cos(theta)*1i)/2-(B.*R(12)*ghx*mub*cos(theta)*1i)/2+(B.*R(20)*ghx*mub*cos(theta)*1i)/2;R(16)*Ua*1i-R(3)*Ub*1i+(R(17)*db*1i)/2-(R(23)*dd*1i)/2-(R(18)*g2)/2-(R(18)*g3)/2-R(18)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+R(18)*(a*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-(B.*R(8)*gex*mub*cos(theta)*1i)/2+(B.*R(20)*gex*mub*cos(theta)*1i)/2-(B.*R(13)*ghx*mub*cos(theta)*1i)/2+(B.*R(19)*ghx*mub*cos(theta)*1i)/2;-R(4)*Ub*1i+R(16)*Ub*1i+(R(20)*dd*1i)/2-(R(24)*dd*1i)/2-R(19)*g3-(B.*R(9)*gex*mub*cos(theta)*1i)/2+(B.*R(17)*gex*mub*cos(theta)*1i)/2-(B.*R(14)*ghx*mub*cos(theta)*1i)/2+(B.*R(18)*ghx*mub*cos(theta)*1i)/2;-R(5)*Ub*1i+(R(19)*dd*1i)/2-(R(25)*dd*1i)/2-(R(20)*g3)/2-(R(20)*g4)/2-R(20)*(a*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+R(20)*(a*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-(B.*R(10)*gex*mub*cos(theta)*1i)/2+(B.*R(18)*gex*mub*cos(theta)*1i)/2-(B.*R(15)*ghx*mub*cos(theta)*1i)/2+(B.*R(17)*ghx*mub*cos(theta)*1i)/2;R(23)*Ua*1i+R(24)*Ub*1i-(R(16)*dd*1i)/2-(R(21)*g4)/2-R(21)*(a*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-(B.*R(11)*gex*mub*cos(theta)*1i)/2-(B.*R(6)*ghx*mub*cos(theta)*1i)/2;(R(23)*db*1i)/2-(R(17)*dd*1i)/2-(R(22)*g1)/2-(R(22)*g4)/2-R(22)*(a*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+R(22)*(a.*B.^2+(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-(B.*R(12)*gex*mub*cos(theta)*1i)/2+(B.*R(24)*gex*mub*cos(theta)*1i)/2-(B.*R(7)*ghx*mub*cos(theta)*1i)/2+(B.*R(25)*ghx*mub*cos(theta)*1i)/2;R(21)*Ua*1i+(R(22)*db*1i)/2-(R(18)*dd*1i)/2-(R(23)*g2)/2-(R(23)*g4)/2-R(23)*(a.*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i+R(23)*(a.*B.^2-(mub*sin(theta)*(gez+ghz).*B)/2+wb-wl)*1i-(B.*R(13)*gex*mub*cos(theta)*1i)/2+(B.*R(25)*gex*mub*cos(theta)*1i)/2-(B.*R(8)*ghx*mub*cos(theta)*1i)/2+(B.*R(24)*ghx*mub*cos(theta)*1i)/2;R(21)*Ub*1i-(R(19)*dd*1i)/2+(R(25)*dd*1i)/2-(R(24)*g3)/2-(R(24)*g4)/2+R(24)*(a.*B.^2-(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-R(24)*(a.*B.^2+(mub*sin(theta)*(gez-ghz).*B)/2+wd-wl)*1i-(B.*R(14)*gex*mub*cos(theta)*1i)/2+(B.*R(22)*gex*mub*cos(theta)*1i)/2-(B.*R(9)*ghx*mub*cos(theta)*1i)/2+(B.*R(23)*ghx*mub*cos(theta)*1i)/2;-(R(20)*dd*1i)/2+(R(24)*dd*1i)/2-R(25)*g4-(B.*R(15)*gex*mub*cos(theta)*1i)/2+(B.*R(23)*gex*mub*cos(theta)*1i)/2-(B.*R(10)*ghx*mub*cos(theta)*1i)/2+(B.*R(22)*ghx*mub*cos(theta)*1i)/2]
tspan = linspace(0, 100, numel(B));
Rm = zeros(numel(tspan), numel(B)); % Preallocate
for k = 1:numel(B)
[t,R]=ode45(@(t,R)dRdt(t,R,B(k)),tspan,[1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]);
Rm(:,k) = R(:,1);
end
ttls = compose('B = %.1f',B);
figure
for k = 1:25
subplot(5,5,k)
% yyaxis left
plot(B,real(Rm(:,1+(k-1)*8)))
% ylabel('Real')
% yyaxis right
% plot(B,imag(Rm(:,1+(k-1)*8)))
% ylabel('Imag')
title(ttls{1+(k-1)*8})
grid
end
Since the plots are a function of the ‘B’ vector, it and ‘tspan’ must have the same number of elements. (If the plots were a function of ‘t’ instead, this would not be necessary.)
.

More Answers (1)

Torsten
Torsten on 15 May 2022
b=[0:0.1:10];
d0=0.2
wb=1
wd=wb-d0
wl=wb
db=0.18
dd=0.05
a=0.02
Ua=0.15
Ub=0.02
g1= 0.001
g2=g1
g3= 0.0001
g4=g3
%syms B theta mub ghx ghz gex gez d0 wb wd wl db dd a Ua Ub g1 g2 g3 g4
theta=pi/3
mub=0.0579; %# Bohr magneton
ghx=-0.35 ;
ghz=-2.2;
gex=-0.65;
gez=-0.8;
for k=1:numel(b)
B = b(k);
bp=mub*B.*sin(theta)*(gez+ghz)*0.5 ; %bp
bm=mub*B.*sin(theta)*(gez-ghz)*0.5; %bm
be=mub*B.*cos(theta)*gex*0.5; %be
bh=mub*B.*cos(theta)*ghx*0.5; %bh
dRdt=@(t,R)[R(7)*g1 + R(13)*g2 + R(19)*g3 + R(25)*g4 + R(3)*Ua*1i - R(11)*Ua*1i + R(4)*Ub*1i - R(16)*Ub*1i
- R(12)*Ua*1i - R(17)*Ub*1i + (R(3)*db*1i)/2 - (R(2)*g1)/2 + R(2)*(a*B.^2 + (mub*sin(theta)*(gez + ghz)*B)/2 + wb - wl)*1i + (B.*R(4)*gex*mub*cos(theta)*1i)/2 + (B.*R(5)*ghx*mub*cos(theta)*1i)/2
R(1)*Ua*1i - R(13)*Ua*1i - R(18)*Ub*1i + (R(2)*db*1i)/2 - (R(3)*g2)/2 + R(3)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(5)*gex*mub*cos(theta)*1i)/2 + (B.*R(4)*ghx*mub*cos(theta)*1i)/2
- R(14)*Ua*1i + R(1)*Ub*1i - R(19)*Ub*1i + (R(5)*dd*1i)/2 - (R(4)*g3)/2 + R(4)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + (B.*R(2)*gex*mub*cos(theta)*1i)/2 + (B.*R(3)*ghx*mub*cos(theta)*1i)/2
- R(15)*Ua*1i - R(20)*Ub*1i + (R(4)*dd*1i)/2 - (R(5)*g4)/2 + R(5)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + (B.*R(3)*gex*mub*cos(theta)*1i)/2 + (B.*R(2)*ghx*mub*cos(theta)*1i)/2
R(8)*Ua*1i + R(9)*Ub*1i - (R(11)*db*1i)/2 - (R(6)*g1)/2 - R(6)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(16)*gex*mub*cos(theta)*1i)/2 - (B.*R(21)*ghx*mub*cos(theta)*1i)/2
(R(8)*db*1i)/2 - (R(12)*db*1i)/2 - R(7)*g1 + (B.*R(9)*gex*mub*cos(theta)*1i)/2 - (B.*R(17)*gex*mub*cos(theta)*1i)/2 + (B.*R(10)*ghx*mub*cos(theta)*1i)/2 - (B.*R(22)*ghx*mub*cos(theta)*1i)/2
R(6)*Ua*1i + (R(7)*db*1i)/2 - (R(13)*db*1i)/2 - (R(8)*g1)/2 - (R(8)*g2)/2 + R(8)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - R(8)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(10)*gex*mub*cos(theta)*1i)/2 - (B*R(18)*gex*mub*cos(theta)*1i)/2 + (B.*R(9)*ghx*mub*cos(theta)*1i)/2 - (B.*R(23)*ghx*mub*cos(theta)*1i)/2
R(6)*Ub*1i - (R(14)*db*1i)/2 + (R(10)*dd*1i)/2 - (R(9)*g1)/2 - (R(9)*g3)/2 + R(9)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(9)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(7)*gex*mub*cos(theta)*1i)/2 - (B.*R(19)*gex*mub*cos(theta)*1i)/2 + (B.*R(8)*ghx*mub*cos(theta)*1i)/2 - (B.*R(24)*ghx*mub*cos(theta)*1i)/2
- (R(15)*db*1i)/2 + (R(9)*dd*1i)/2 - (R(10)*g1)/2 - (R(10)*g4)/2 + R(10)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(10)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(8)*gex*mub*cos(theta)*1i)/2 - (B.*R(20)*gex*mub*cos(theta)*1i)/2 + (B.*R(7)*ghx*mub*cos(theta)*1i)/2 - (B.*R(25)*ghx*mub*cos(theta)*1i)/2
- R(1)*Ua*1i + R(13)*Ua*1i + R(14)*Ub*1i - (R(6)*db*1i)/2 - (R(11)*g2)/2 - R(11)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(21)*gex*mub*cos(theta)*1i)/2 - (B.*R(16)*ghx*mub*cos(theta)*1i)/2
- R(2)*Ua*1i - (R(7)*db*1i)/2 + (R(13)*db*1i)/2 - (R(12)*g1)/2 - (R(12)*g2)/2 - R(12)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + R(12)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(14)*gex*mub*cos(theta)*1i)/2 - (B.*R(22)*gex*mub*cos(theta)*1i)/2 + (B.*R(15)*ghx*mub*cos(theta)*1i)/2 - (B.*R(17)*ghx*mub*cos(theta)*1i)/2
- R(3)*Ua*1i + R(11)*Ua*1i - (R(8)*db*1i)/2 + (R(12)*db*1i)/2 - R(13)*g2 + (B.*R(15)*gex*mub*cos(theta)*1i)/2 - (B.*R(23)*gex*mub*cos(theta)*1i)/2 + (B.*R(14)*ghx*mub*cos(theta)*1i)/2 - (B.*R(18)*ghx*mub*cos(theta)*1i)/2
- R(4)*Ua*1i + R(11)*Ub*1i - (R(9)*db*1i)/2 + (R(15)*dd*1i)/2 - (R(14)*g2)/2 - (R(14)*g3)/2 + R(14)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(14)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(12)*gex*mub*cos(theta)*1i)/2 - (B.*R(24)*gex*mub*cos(theta)*1i)/2 + (B.*R(13)*ghx*mub*cos(theta)*1i)/2 - (B.*R(19)*ghx*mub*cos(theta)*1i)/2
- R(5)*Ua*1i - (R(10)*db*1i)/2 + (R(14)*dd*1i)/2 - (R(15)*g2)/2 - (R(15)*g4)/2 + R(15)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(15)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i + (B.*R(13)*gex*mub*cos(theta)*1i)/2 - (B.*R(25)*gex*mub*cos(theta)*1i)/2 + (B.*R(12)*ghx*mub*cos(theta)*1i)/2 - (B.*R(20)*ghx*mub*cos(theta)*1i)/2
R(18)*Ua*1i - R(1)*Ub*1i + R(19)*Ub*1i - (R(21)*dd*1i)/2 - (R(16)*g3)/2 - R(16)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(6)*gex*mub*cos(theta)*1i)/2 - (B.*R(11)*ghx*mub*cos(theta)*1i)/2
- R(2)*Ub*1i + (R(18)*db*1i)/2 - (R(22)*dd*1i)/2 - (R(17)*g1)/2 - (R(17)*g3)/2 - R(17)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(17)*(a*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(7)*gex*mub*cos(theta)*1i)/2 + (B.*R(19)*gex*mub*cos(theta)*1i)/2 - (B.*R(12)*ghx*mub*cos(theta)*1i)/2 + (B.*R(20)*ghx*mub*cos(theta)*1i)/2
R(16)*Ua*1i - R(3)*Ub*1i + (R(17)*db*1i)/2 - (R(23)*dd*1i)/2 - (R(18)*g2)/2 - (R(18)*g3)/2 - R(18)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(18)*(a*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(8)*gex*mub*cos(theta)*1i)/2 + (B.*R(20)*gex*mub*cos(theta)*1i)/2 - (B.*R(13)*ghx*mub*cos(theta)*1i)/2 + (B.*R(19)*ghx*mub*cos(theta)*1i)/2
- R(4)*Ub*1i + R(16)*Ub*1i + (R(20)*dd*1i)/2 - (R(24)*dd*1i)/2 - R(19)*g3 - (B.*R(9)*gex*mub*cos(theta)*1i)/2 + (B.*R(17)*gex*mub*cos(theta)*1i)/2 - (B.*R(14)*ghx*mub*cos(theta)*1i)/2 + (B.*R(18)*ghx*mub*cos(theta)*1i)/2
- R(5)*Ub*1i + (R(19)*dd*1i)/2 - (R(25)*dd*1i)/2 - (R(20)*g3)/2 - (R(20)*g4)/2 - R(20)*(a*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(20)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(10)*gex*mub*cos(theta)*1i)/2 + (B.*R(18)*gex*mub*cos(theta)*1i)/2 - (B.*R(15)*ghx*mub*cos(theta)*1i)/2 + (B.*R(17)*ghx*mub*cos(theta)*1i)/2
R(23)*Ua*1i + R(24)*Ub*1i - (R(16)*dd*1i)/2 - (R(21)*g4)/2 - R(21)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(11)*gex*mub*cos(theta)*1i)/2 - (B.*R(6)*ghx*mub*cos(theta)*1i)/2
(R(23)*db*1i)/2 - (R(17)*dd*1i)/2 - (R(22)*g1)/2 - (R(22)*g4)/2 - R(22)*(a*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(22)*(a.*B.^2 + (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(12)*gex*mub*cos(theta)*1i)/2 + (B.*R(24)*gex*mub*cos(theta)*1i)/2 - (B.*R(7)*ghx*mub*cos(theta)*1i)/2 + (B.*R(25)*ghx*mub*cos(theta)*1i)/2
R(21)*Ua*1i + (R(22)*db*1i)/2 - (R(18)*dd*1i)/2 - (R(23)*g2)/2 - (R(23)*g4)/2 - R(23)*(a.*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i + R(23)*(a.*B.^2 - (mub*sin(theta)*(gez + ghz).*B)/2 + wb - wl)*1i - (B.*R(13)*gex*mub*cos(theta)*1i)/2 + (B.*R(25)*gex*mub*cos(theta)*1i)/2 - (B.*R(8)*ghx*mub*cos(theta)*1i)/2 + (B.*R(24)*ghx*mub*cos(theta)*1i)/2
R(21)*Ub*1i - (R(19)*dd*1i)/2 + (R(25)*dd*1i)/2 - (R(24)*g3)/2 - (R(24)*g4)/2 + R(24)*(a.*B.^2 - (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - R(24)*(a.*B.^2 + (mub*sin(theta)*(gez - ghz).*B)/2 + wd - wl)*1i - (B.*R(14)*gex*mub*cos(theta)*1i)/2 + (B.*R(22)*gex*mub*cos(theta)*1i)/2 - (B.*R(9)*ghx*mub*cos(theta)*1i)/2 + (B.*R(23)*ghx*mub*cos(theta)*1i)/2
- (R(20)*dd*1i)/2 + (R(24)*dd*1i)/2 - R(25)*g4 - (B.*R(15)*gex*mub*cos(theta)*1i)/2 + (B.*R(23)*gex*mub*cos(theta)*1i)/2 - (B.*R(10)*ghx*mub*cos(theta)*1i)/2 + (B.*R(22)*ghx*mub*cos(theta)*1i)/2]
[t{k},R{k}]=ode45(dRdt,[0 100],[1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0])
end
  2 Comments
Celso Júnior
Celso Júnior on 22 May 2022
Good Morning!
I'm trying to plot the solutions of the differential equations of this attached system, as a function of the magnetic field, in the same way you did in the previous calculation. I wanted to plot plot(t,R:,19) with several different values ​​of B B=0:1:10. You are helping me a lot and I already appreciate your help and patience.
B=7
d0=0.2
wb=1
wd=wb-d0
wl=wb
db=0.18
dd=0.05
a=0.02
Ua=1.5 %0.15
Ub=2 %0.02
g1= 0.1 %0.001
g2=g1
g3= 0.1 %0.0001
g4=g3
%syms B theta mub ghx ghz gex gez d0 wb wd wl db dd a Ua Ub g1 g2 g3 g4
theta=pi/3
mub=0.0579; %# Bohr magneton
ghx=-0.35 ;
ghz=-2.2;
gex=-0.65;
gez=-0.8;
bp=mub*B.*sin(theta)*(gez+ghz)*0.5 ; %bp
bm=mub*B.*sin(theta)*(gez-ghz)*0.5; %bm
be=mub*B.*cos(theta)*gex*0.5; %be
bh=mub*B.*cos(theta)*ghx*0.5; %bh
dRdt=@(t,R)[ R(7)*g1 + R(13)*g2 + R(19)*g3 + R(25)*g4 + R(2)*Ua*1i - R(6)*Ua*1i + R(3)*Ub*1i - R(11)*Ub*1i
R(1)*Ua*1i - R(7)*Ua*1i - R(12)*Ub*1i + R(4)*be*1i + R(5)*bh*1i + (R(3)*db*1i)/2 - (R(2)*g1)/2 + R(2)*(a*B.^2 + bp + wb - wl)*1i
- R(8)*Ua*1i + R(1)*Ub*1i - R(13)*Ub*1i + R(5)*be*1i + R(4)*bh*1i + (R(2)*db*1i)/2 - (R(3)*g2)/2 - R(3)*(- a*B.^2 + bp - wb + wl)*1i
- R(9)*Ua*1i - R(14)*Ub*1i + R(2)*be*1i + R(3)*bh*1i + (R(5)*dd*1i)/2 - (R(4)*g3)/2 - R(4)*(- a*B.^2 + bm - wd + wl)*1i
- R(10)*Ua*1i - R(15)*Ub*1i + R(3)*be*1i + R(2)*bh*1i + (R(4)*dd*1i)/2 - (R(5)*g4)/2 + R(5)*(a*B.^2 + bm + wd - wl)*1i
- R(1)*Ua*1i + R(7)*Ua*1i + R(8)*Ub*1i - R(16)*be*1i - R(21)*bh*1i - (R(11)*db*1i)/2 - (R(6)*g1)/2 - R(6)*(a*B.^2 + bp + wb - wl)*1i
- R(2)*Ua*1i + R(6)*Ua*1i + R(9)*be*1i - R(17)*be*1i + R(10)*bh*1i - R(22)*bh*1i + (R(8)*db*1i)/2 - (R(12)*db*1i)/2 - R(7)*g1
- R(3)*Ua*1i + R(6)*Ub*1i + R(10)*be*1i - R(18)*be*1i + R(9)*bh*1i - R(23)*bh*1i + (R(7)*db*1i)/2 - (R(13)*db*1i)/2 - (R(8)*g1)/2 - (R(8)*g2)/2 - R(8)*(a*B.^2 + bp + wb - wl)*1i - R(8)*(- a*B.^2 + bp - wb + wl)*1i
- R(4)*Ua*1i + R(7)*be*1i - R(19)*be*1i + R(8)*bh*1i - R(24)*bh*1i - (R(14)*db*1i)/2 + (R(10)*dd*1i)/2 - (R(9)*g1)/2 - (R(9)*g3)/2 - R(9)*(- a*B.^2 + bm - wd + wl)*1i - R(9)*(a*B.^2 + bp + wb - wl)*1i
- R(5)*Ua*1i + R(8)*be*1i - R(20)*be*1i + R(7)*bh*1i - R(25)*bh*1i - (R(15)*db*1i)/2 + (R(9)*dd*1i)/2 - (R(10)*g1)/2 - (R(10)*g4)/2 + R(10)*(a*B.^2 + bm + wd - wl)*1i - R(10)*(a*B.^2 + bp + wb - wl)*1i
R(12)*Ua*1i - R(1)*Ub*1i + R(13)*Ub*1i - R(21)*be*1i - R(16)*bh*1i - (R(6)*db*1i)/2 - (R(11)*g2)/2 + R(11)*(- a*B.^2 + bp - wb + wl)*1i
R(11)*Ua*1i - R(2)*Ub*1i + R(14)*be*1i - R(22)*be*1i + R(15)*bh*1i - R(17)*bh*1i - (R(7)*db*1i)/2 + (R(13)*db*1i)/2 - (R(12)*g1)/2 - (R(12)*g2)/2 + R(12)*(a*B.^2 + bp + wb - wl)*1i + R(12)*(- a*B.^2 + bp - wb + wl)*1i
- R(3)*Ub*1i + R(11)*Ub*1i + R(15)*be*1i - R(23)*be*1i + R(14)*bh*1i - R(18)*bh*1i - (R(8)*db*1i)/2 + (R(12)*db*1i)/2 - R(13)*g2
- R(4)*Ub*1i + R(12)*be*1i - R(24)*be*1i + R(13)*bh*1i - R(19)*bh*1i - (R(9)*db*1i)/2 + (R(15)*dd*1i)/2 - (R(14)*g2)/2 - (R(14)*g3)/2 - R(14)*(- a*B.^2 + bm - wd + wl)*1i + R(14)*(- a*B.^2 + bp - wb + wl)*1i
- R(5)*Ub*1i + R(13)*be*1i - R(25)*be*1i + R(12)*bh*1i - R(20)*bh*1i - (R(10)*db*1i)/2 + (R(14)*dd*1i)/2 - (R(15)*g2)/2 - (R(15)*g4)/2 + R(15)*(a*B.^2 + bm + wd - wl)*1i + R(15)*(- a*B.^2 + bp - wb + wl)*1i
R(17)*Ua*1i + R(18)*Ub*1i - R(6)*be*1i - R(11)*bh*1i - (R(21)*dd*1i)/2 - (R(16)*g3)/2 + R(16)*(- a*B.^2 + bm - wd + wl)*1i
R(16)*Ua*1i - R(7)*be*1i + R(19)*be*1i - R(12)*bh*1i + R(20)*bh*1i + (R(18)*db*1i)/2 - (R(22)*dd*1i)/2 - (R(17)*g1)/2 - (R(17)*g3)/2 + R(17)*(- a*B.^2 + bm - wd + wl)*1i + R(17)*(a*B.^2 + bp + wb - wl)*1i
R(16)*Ub*1i - R(8)*be*1i + R(20)*be*1i - R(13)*bh*1i + R(19)*bh*1i + (R(17)*db*1i)/2 - (R(23)*dd*1i)/2 - (R(18)*g2)/2 - (R(18)*g3)/2 + R(18)*(- a*B.^2 + bm - wd + wl)*1i - R(18)*(- a*B.^2 + bp - wb + wl)*1i
- R(9)*be*1i + R(17)*be*1i - R(14)*bh*1i + R(18)*bh*1i + (R(20)*dd*1i)/2 - (R(24)*dd*1i)/2 - R(19)*g3
- R(10)*be*1i + R(18)*be*1i - R(15)*bh*1i + R(17)*bh*1i + (R(19)*dd*1i)/2 - (R(25)*dd*1i)/2 - (R(20)*g3)/2 - (R(20)*g4)/2 + R(20)*(a*B.^2 + bm + wd - wl)*1i + R(20)*(- a*B.^2 + bm - wd + wl)*1i
R(22)*Ua*1i + R(23)*Ub*1i - R(11)*be*1i - R(6)*bh*1i - (R(16)*dd*1i)/2 - (R(21)*g4)/2 - R(21)*(a*B.^2 + bm + wd - wl)*1i
R(21)*Ua*1i - R(12)*be*1i + R(24)*be*1i - R(7)*bh*1i + R(25)*bh*1i + (R(23)*db*1i)/2 - (R(17)*dd*1i)/2 - (R(22)*g1)/2 - (R(22)*g4)/2 - R(22)*(a*B.^2 + bm + wd - wl)*1i + R(22)*(a*B.^2 + bp + wb - wl)*1i
R(21)*Ub*1i - R(13)*be*1i + R(25)*be*1i - R(8)*bh*1i + R(24)*bh*1i + (R(22)*db*1i)/2 - (R(18)*dd*1i)/2 - (R(23)*g2)/2 - (R(23)*g4)/2 - R(23)*(a*B.^2 + bm + wd - wl)*1i - R(23)*(- a*B.^2 + bp - wb + wl)*1i
- R(14)*be*1i + R(22)*be*1i - R(9)*bh*1i + R(23)*bh*1i - (R(19)*dd*1i)/2 + (R(25)*dd*1i)/2 - (R(24)*g3)/2 - (R(24)*g4)/2 - R(24)*(a*B.^2 + bm + wd - wl)*1i - R(24)*(- a*B.^2 + bm - wd + wl)*1i
- R(15)*be*1i + R(23)*be*1i - R(10)*bh*1i + R(22)*bh*1i - (R(20)*dd*1i)/2 + (R(24)*dd*1i)/2 - R(25)*g4]
[t,R]=ode45(dRdt,[0 150],[1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]);
plot(t,R(:,19))
legend('\rho_{33}')
Torsten
Torsten on 22 May 2022
Edited: Torsten on 22 May 2022
Did you try the above solution ? I don't see that your question differs from the first.
Since the solution is complex, plotting in your case means plotting real(R ), imag(R), abs(R ), angle(R ) ...

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!