Editing Talk:Forum for GNU Octave

Dear All,
I am trying to write an transient analysis solver for a multi degree of freedom system. it consists of mass, stiffness, damping, gyroscopic damping matrices and so on. my initial condition is declared by zeros(2*nr,1). after I run the code, the output y which is a matrix of 6x41  has all it's rows and columns just zero. It would be really apprectiated if some one has an advice / opinion about what am I doing wrong?
for the main .m file it is as below:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = [0.0020*2*pi 8*pi 0];
tspan = [0 4] ;
frunb=zeros(size(K,1),1); %%% unbalance vector
frunb(5,1) = 1; %%% 
nr = 3;
options = odeset;
[t,y] = ode45(@deriv,tspan,zeros(2*nr,1),options,M,K,C,G,alpha,frunb);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


and below is the function deriv (saved in a different file named as deriv.m)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dz = deriv(t,z,M,K,C,G,alpha,ubforce)
  
  phi = alpha(1)*t*t + alpha(2)*t + alpha(3);
  d_phi = 2*alpha(1)*t + alpha(2);
  dd_phi = 2*alpha(1);
  A = [(-inv(M)*(C-j*phi*G)) (-inv(M)*K); eye(size(M)) zeros(size(M))];
  [U,D] = eig(A);
  uncoupled_A = (inv(U))*A*(U);
  B = [inv(M)*ubforce;zeros(size(M))*ubforce];
  uncoupled_B = inv(U) * B;
  double_uncoupled_A_1 = [(uncoupled_A(1,1) + uncoupled_A(2,2)) (-uncoupled_A(1,1)*uncoupled_A(2,2)); 1 0]; 
  double_uncoupled_A_2 = [(uncoupled_A(3,3) + uncoupled_A(4,4)) (-uncoupled_A(3,3)*uncoupled_A(4,4)); 1 0];
  double_uncoupled_A_3 = [(uncoupled_A(5,5) + uncoupled_A(6,6)) (-uncoupled_A(5,5)*uncoupled_A(6,6)); 1 0];
  U11 = [uncoupled_A(1,1)  uncoupled_A(2,2); 1 1];
  U22 = [uncoupled_A(3,3)  uncoupled_A(4,4); 1 1];
  U33 = [uncoupled_A(5,5)  uncoupled_A(6,6); 1 1];
  U = U(1:6,1:6);
  U1 = blkdiag(U11,U22,U33);
  double_uncoupled_A = blkdiag(double_uncoupled_A_1,double_uncoupled_A_2,double_uncoupled_A_3);
  double_uncoupled_B = U1 * inv(U) * B(1:6,:);
  gain = [1 0 0 0 0 0] * (phi^2)*exp(j*phi*t)*z;
  dz = (double_uncoupled_A * z) + double_uncoupled_B * gain';
endfunction