ode45 funtion部分中微分表达式中的参数仍是时变的怎么办

clc
clear all

a=-1;
b=0;
c=b-a;

I11=200;
I22=150;
I33=120;
T1=2;
T2=1.25;
T3=5;
P1=2I11/T1;
P2=2I22/T2;
P3=2I33/T3;
K1=P1^2(1-a)^2/2/I11;
K2=P1^2*(1-a)^2/2/I11;
K3=P1^2*(1-a)^2/2/I11;

tspan=[0 200];
y0=[1;1;1;0;0;0;0]; %y=[w(1);w(2);w(3);fai;deta(1);deta(2);deta(3)];
[t,y]=ode45(@hanshu1,tspan,y0);

deta=sqrt(deta1^2+deta2^2+deta3^2);
sum1=adeta^2+sqrt(c^2+cdeta^2*(1-a^2));
sum2=(-ac+sqrt(c^2+cdeta^2*(1-a^2)))/c;
w=sqrt(w1^2+w2^2+w3^2);
u1=(I33-I22)w2w3-K1*(1+a/sum2)deta1-P1w1;
u2=(I11-I33)w3w1-K2*(1+a/sum2)deta2-P2w2;
u3=(I22-I11)w1w2-K3*(1+a/sum2)deta3-P3w3;

figure(1)
plot(t,y(:,1),'r',t,y(:,2),'g',t,y(:,3),'b');
legend('w1','w2','w3');
title('角速度');
xlabel('t/s');
ylabel('rad/s');
hold on

figure(2)
plot(t,y(:,4));
title('主旋转角度');
xlabel('t/s');
ylabel('deg');
hold on

figure(3)
plot(t,y(:,5),'r',t,y(:,6),'g',t,y(:,7),'b');
legend('deta1','deta2','deta3');
title('姿态');
xlabel('t/s');
hold on

figure(4)
plot(t,u1,'r',t,u2,'g',t,u3,'b');
legend('u1','u2','u3');
title('控制力矩');
xlabel('t/s');
ylabel('N*M');
grid on

function dy = hanshu1(t,y)

w1=y(1);
w2=y(2);
w3=y(3);

fai=y(4);

deta1=y(5);
deta2=y(6);
deta3=y(7);

dw(1)=(-(I33-I22)w2w3)/I11+u1;
dw(2)=(-(I11-I33)w3w1)/I22+u2;
dw(3)=(-(I22-I11)w1w2)/I33+u3;

dfai=2*w;

ddeta(1)=((sum1/sum2+deta1^2)w1+(deta1deta2-deta3)w2+(deta1deta3+deta2)w3)/2;
ddeta(2)=((deta1deta2+deta3)*w1+(sum1/sum2+deta2^2)w2+(deta2deta3-deta1)w3)/2;
ddeta(3)=((deta1deta3-deta2)w1+(deta2deta3+deta1)*w2+(sum1/sum2+deta3^2)*w3)/2;

dy=[dw(1);dw(2);dw(3);dfai;ddeta(1);ddeta(2);ddeta(3)];

end

其中deta=sqrt(deta1^2+deta2^2+deta3^2);
sum1=adeta^2+sqrt(c^2+cdeta^2*(1-a^2));
sum2=(-ac+sqrt(c^2+cdeta^2*(1-a^2)))/c;
w=sqrt(w1^2+w2^2+w3^2);
u1=(I33-I22)w2w3-K1*(1+a/sum2)deta1-P1w1;
u2=(I11-I33)w3w1-K2*(1+a/sum2)deta2-P2w2;
u3=(I22-I11)w1w2-K3*(1+a/sum2)deta3-P3w3;
这一部分不知道怎么处理 因为ode45积分的w1，w2，w3，fai，deta1，deta2，deta3表达式中包括w，deta，sum1，sum2这些也是时变 是应该嵌套还是再次ode45 求解释

把w1，w2，w3，fai，deta1，deta2，deta3关于w，deta，sum1，sum2的表达式写在hanshu1函数中，这些状态都可以是时变的，ode45能够求解的。