Error，Need some help

In order to solve the equations， I edited 2 M files. In equation (6), X is coefficient which varies with t(t is time)

【ODE451_main】

clc;clear

tspan=[0,180];

y0=[1e-5;0;0;0;0];

[t,y]=ode23('ODE45_fun',tspan,y0);

[m,n]=size(y);

for i=1:m

      t=[0 25 50 70 80 100 180];

x=[0 5  7  9  10 11 14];

%plot(t,x,'r*'); hold on;

k = 3;  %多项式的阶次

p=polyfit(t,x,k);   %得到拟合的多项式系数

%dt = linspace(0,100,100);
%dx = polyval(p,dt);
%plot(dt,dx);
%xlabel('t');
%ylabel('x');
%title('x-t的拟合曲线');

syms t;%syms t Y;

X = [t^3 t^2 t 1]* p'; %得到X-t的函数

z =-178075801.6*y(i,4)*y(i,5)+X*y(i,3)*y(i,2)+1119396.815471592+7481.104004*sin(2.047*t(i)); %得到z(t,Y）的符号表达式

Z = matlabFunction(z);   %将符号表达式转化为函数

t = linspace(0,180,180);

y(i,6) = Z(t);  %Dz-t的方程;dz = Z(3,t);y(i,6) = Z(t); 

end
data=[t,y]; 

save ODE45_data.txt data -ascii

subplot(2,3,1),plot(t,y(:,1)),title('y(1)')

xlabel('t');ylabel('y');

subplot(2,3,2),plot(t,y(:,2)),title('y(2)')

xlabel('t');ylabel('y');

subplot(2,3,3),plot(t,y(:,3)),title('y(3)')

xlabel('t');ylabel('y');

subplot(2,3,4),plot(t,y(:,4)),title('y(4)')

xlabel('t');ylabel('y');

subplot(2,3,5),plot(t,y(:,5)),title('y(5)')

xlabel('t');ylabel('y');

subplot(2,3,6),plot(t,y(:,6)),title('y(6)')

xlabel('t');ylabel('y');

% plot(t,y(:,1),'bo',t,y(:,2),'rx',t,y(:,3),'gv',t,y(:,4),'r-');
grid on

【ODE45_fun】

function  dy=ODE45_fun(t,y)

dy(1)=-1.918298553*y(3)*y(4)-121.6697369*y(5)*y(2)+0.006472085*y(2)*y(2)+15.25250926*y(5)*y(5)-0.518363603*sin(2.047*t)+0.001124759;

dy(2)=0.007229182*y(5)*y(1)-0.013867729*y(3)-0.005151943*sin(2.047*t)+33.43424564*y(4)*y(5)-0.092169794*y(3)*y(2)-0.698266828;

dy(3)=72.10986245*y(2)+0.52129529*y(4)*y(1)+0.025471074*y(3)+921.886526*y(4)/y(1)-0.47870471*y(4)+0.025471074*y(1)*t-0.38220722*cos(2.047*t)-4.62279911;

dy(4)=-57.37263009*y(5)+0.001053501*y(3)*y(1)+20.28825016/y(1)-0.001053501*y(3)+0.064741605*y(4)+20.28825016*t-0.006201915*sin(2.047*t)-3651.885374303154;

dy(5)=0.017284293*y(4)+0.00278644*y(1)*y(2)-0.551218454*y(2)*y(5)*y(5)+0.010839281*y(2)*y(2)*y(5)+0.020353114*cos(2.047*t)+0.110594984;

dy=[dy(1);dy(2);dy(3);dy(4);dy(5)];

[Walter Roberson]

Change your line

y(i,6) = Z(t);

to

Zt = Z(t);
fprintf('iteration i = %d, number of elements in Zt is %d\n', i, numel(Zt));
y(i,6) = Zt;

This will not prevent the error but it will help you track down the circumstances under which the error occurs. Does it occur on the first iteration, or does it occur after it has been running for some time?