how can i use ode45 to solve a matrix form differential equation?

hi every one there are 4 differential eqs that i want to use ode45 to solve them,but one of eqs is a 3*1 column matrix eq. because of it,i encounter with this error: ??? Index exceeds matrix dimensions.

Error in ==> zahra_thesis921_attitude_hold_simulation_ode at 36 wtilda=[0 -w(3,:) w(2,:);w(3,:) 0 -w(1,:);-w(2,:) w(1,:) 0];

I have never used ode45 for a matrix form eq,so i guess this is the reason of error. please give me some advices to appropriately solve this problem and correctly run it. thanks alot... here are my codes: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the first code is constants and initial conditions that named: "zahra_thesis921_attitude_hold_simulation_input": % zahra_thesis921_attitude_hold_simulation_input

clc; close all clear all global psi0 theta0 phi0 w0 t0 tf %%%%%% Input File Call

psi0=2;theta0=2;phi0=2; %deg w0=[0.2;0.2;0.2]; %deg/sec

t0=0; tf=50; % final time %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the second code is my ode eq named "zahra_thesis921_attitude_hold_simulation_ode":

function ydot=zahra_thesis921_attitude_hold_simulation_ode(t,y) global psi0 theta0 phi0 w0 t0 tf

w=y(1); psi=y(2); theta=y(3); phi=y(4);

wtilda=[0 -w(3,:) w(2,:);w(3,:) 0 -w(1,:);-w(2,:) w(1,:) 0]; I=diag([39.6,55,55]); S=zeros(6,3); E=eye(6); A(1:6,1:3)=zeros(6,3);A(1:3,4:6)=eye(3);A(4:6,4:6)=zeros(3,3); B(1:3,1:3)=zeros(3,3);B(4:6,1:3)=I^-1; Q(1:3,1:3)=eye(3);Q(1:3,4:6)=zeros(3,3);Q(4:6,1:3)=zeros(3,3);Q(4:6,4:6)=zeros(3,3); R=(2*10^-5)*eye(3); [M,L,G] = care(A,B,Q,R,S,E); Kss=inv®*transp(B)*M x=[psi;theta;phi;w(3,:);w(2,:);w(1,:)]; u=-Kss*x;

% ODEs:

wdot=-(I^-1)*wtilda*I*w+(I^-1)*u; %kinematic differential equations for the chosen 3-2-1 Euler angle set are: psidot=(w(2,:)*sin(phi)+w(3,:)*cos(phi))/cos(theta); thetadot=w(2,:)*cos(phi)-w(3,:)*sin(phi); phidot=w(1,:)+(w(2,:)*sin(phi)+w(3,:)*cos(phi))*tan(theta);

ydot=[wdot psidot thetadot phidot]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the third and last code is "zahra_thesis921_attitude_hold_simulation_main" that you must run it: global psi0 theta0 phi0 w0 t0 tf

zahra_thesis921_attitude_hold_simulation_input;

tspan=[t0 tf]; % X0: ic=[w0 num2cell(psi0) num2cell(theta0) num2cell(phi0)];

[t,Y]=ode45(@zahra_thesis921_attitude_hold_simulation_ode,tspan,ic); w=Y(:,1); %size(Y); psi=Y(:,2); theta=Y(:,3); phi=Y(:,4); %% control effort x=[psi;theta;phi;w(3,:);w(2,:);w(1,:)]; Kss=[223.6068,0,0,133.0776,0,0;0,223.6068,0,0,156.8335,0;0,0,223.6068,0,0,156.8335]; u=-Kss.*x;

figure(1) % Plot angular velocity subplot(3,1,1);plot(t,w(1,:));xlabel('t');ylabel('wx<rad/s>') grid; subplot(3,1,2);plot(t,w(2,:));xlabel('t');ylabel('wy<rad/s>') grid; subplot(3,1,3);plot(t,w(3,:));xlabel('Time sec');ylabel('wz<rad/s>') grid; figure(2) % Plot Euler angle subplot(3,1,1);plot(t,psi);xlabel('t');ylabel('\psi<rad>') grid; subplot(3,1,2);plot(t,theta);xlabel('t');ylabel('\theta<rad>') grid; subplot(3,1,3);plot(t,phi);xlabel('Time sec');ylabel('\phi<rad>') grid;

% control effort figure(3) plot(t,u(1,:),t,u(2,:),'--',t,u(3,:),'-.');xlabel('t');ylabel('control Input/Torque-foot-lbs') title('Control Torque History')