[t,y]=ode45('rigid',[1 10],[2 -3]);
 
figure(1);
clf;
plot(t,y(:,1),'r-');
hold on;
z=1./t+1./t.^2;
plot(t,z,'k:');
