function dy = rigid(t,y)
     dy = zeros(2,1);    % a column vector
     dy(1) = y(2);
     dy(2) = -2*y(1)/(t^2)-4*y(2)/t;
