// u'' + 0.125 u' + 4 u = 3 cos ( omega t)
// u' = z
function wdot = f(t, w)
wdot = [ w(2,:); -0.125*w(2,:) - 4*w(1,:)+3*cos(omega*t) ];
endfunction
omega=2; t0=0; w0=[0;2];t=0:0.1:40*%pi;
w=ode(w0,t0,t,f);
scf(0);
clf()
plot(t,3*cos(omega*t),"g",t,w(1,:),"r");
xtitle("input (green) equal omega_0 -- output red","time t","u");
scf(1);
clf()
plot(t,3*cos(omega*t),"g",t,w(1,:),"r");
xtitle("input (green) equal omega_0 -- output red","time t","u");
zoom_rect([30*%pi,-15,35*%pi,15]);