function harmonic;
global w2 b wf af;
w2=4.1
wf=2.
b=0.1
af=0.2
x0=[0.;0.];
[t,x]=ode45(@rpp,[0 100], x0);
plot (x(:,1))
function f=rpp(t,x)
global w2 b wf af;
f=[x(2);-w2*x(1)-b*x(2)-af*cos(wf*t)]
