%%% odey1_case.m %%%%
global e xi Omega omega
e=0.1;xi=0.01;Omega=0.5;omega=1; t=0:0.1:10;
[t,y]=ode45('odey1',t,[0 0]); plot(t,y(:,1),'r',t,y(:,2),'b')
function dy=odey1(t,y)
%%% odey1.m %%%
global e xi Omega omega
dy=[y(2);e*sin(omega*t)-2*xi/Omega*y(2)-1/Omega^2*y(1)];
end