%
%         spring_rhs
%
  function dy = spring_rhs(t,y,params);

  K = params(1);
  C = params(2);

  z = y(1);
  z_dot = y(2);
  z_k = y(3);
  z_k_dot = y(4);
  z_c = y(5);
  z_c_dot = y(6);

  z_vec = [z; z_dot; z_k; z_k_dot; z_c; z_c_dot];
  J = [0 1; -K -C];
  Der = blkdiag(J,J,J);
  rhs = [0; 0; 0; -z; 0; -z_dot];

%
%
  dy = Der*z_vec+ rhs;