Since your matrices are of 3x3 size the result should be of 3x1 size
function dydt = dynamics(t,y,Mt,Ct,Kt)
dydt = zeros(6,1);
y = y(:); % make column;
dydt(1:3) = y(4:6); % xd = [hd; thetad; betad]
dydt(4:6) = inv(Mt)*(-real(Ct)*y(4:6)-real(Kt)*y(1:3)); % xdd = [hdd; thetadd;betadd]
end
You have 3 variable of second degree each. YOu should have 6 initial conditions
X0 = [h;theta;beta;hd;thetad;betad];
[t,Y] = ode45(@(t,Y)dynamics(t,Y,Mt,Ct,Kt),t,X0);
