function output=thefunc(t,y) % yp(1)=y(2)-y(3)+t; % yp(2)=3*t^2; % yp(3)=y(2)+exp(-t); beta=1; alpha=.005; r=.5; yp(1)=y(1); yp(2)=beta*cos(r*t)-alpha*y(2)-y(1)-y(1)^3; output=yp'; end