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)
Show older comments
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))
0 Comments
Accepted Answer
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
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,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.)
.
2 Comments
More Answers (1)
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
See Also
Categories
Find more on Transformers in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!