How to build a multigraph by varying a variable?

Good evening! I'm trying to build a multigraph, where each graph would have a different theta value. I'm varying theta from 0 to pi . Follow code. Thank you very much in advance.

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 %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=0:0.1:pi

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,theta)[ 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 6000],[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.5;0;0;0;0;0;0.5]);

plot(t,R(:,1),t,R(:,7),t,R(:,13),t,R(:,19),t,R(:,25))

legend('\rho_{00}','\rho_{11}','\rho_{22}','\rho_{33}','\rho_{44}')