The subscripted assignment dimension error is caused when you try to assign a vector of length 1651x1 to a 1x1 subarray of z_dot. The problem is that ode45 does not know that psi_dot and psi_ddot are to be interpreted as a time series, so it computes
c1.*pd.*pd*sin(2*z(1))...
+c2.*pdd*cos(z(1))...
+c3.*sign(pd).*pd.^2*(3.48*cos(z(1))).*(0.82*abs(sign(pd)*z(1)+pi/2)/pi+0.05)...
+c4*abs(z(2))*z(2)...
+c5*z(1);
where pd and pdd are the full 1651x1 vectors psi_dot and psi_ddot.
To correct this, you need to compute pd and pdd from the data and from t. One way is as follows:
First, edit the definition of MYODE to accept an additional input tpsi which is a vector of times corresponding to the elements in psi_dot and psi_ddot. Then interpolate (I have used linear interpolation, but you could use nearest-neighbor interpolation or something else as appropriate) to get pd and pdd:
function [z_dot] = MYODE(t,z,c,tpsi,psi_dot,psi_ddot)
% z(1) = theta % z(2) = thetadot
% time-varying data
pd = interp1(tpsi, psi_dot, t, 'linear');
pdd = interp1(tpsi, psi_ddot, t, 'linear');
% Constants:
c1 = c(1);
c2 = c(2);
c3 = c(3);
c4 = c(4);
c5 = c(5);
% thetadot
z_dot(1,1) = z(2);
% thetadot_dot
z_dot(2,1) = c1.*pd.*pd*sin(2*z(1))...
+c2.*pdd*cos(z(1))...
+c3.*sign(pd).*pd.^2*(3.48*cos(z(1))).*(0.82*abs(sign(pd)*z(1)+pi/2)/pi+0.05)...
+c4*abs(z(2))*z(2)...
+c5*z(1);
end
Then, in your main script or function, you just need to pass in tpsi:
tpsi = (0:1/6600:0.25)';
step = 1650; % N/A
t0 = 0; % Start time (sec)
t1 = 0.25; % Final time (sec)
t = linspace(t0 , t1 , (step+1));
[t,z1] = ode45(@(t,z) MYODE(t,z,c0,tpsi, psi_dot, psi_ddot),t,z0);
