function y=f(t,x)
 y=x*(1-x/(1+eps*cos(om*t)))
endfunction

res=x_dialog("dare om e eps",["5.000","0.100"]);
om=evstr(res(1)); eps=evstr(res(2)); t0=0;

function y=p(x)
  n=length(x)
  for i=1:n
    xin=x(i)
    z=ode(xin,t0,[t0, t0+2*%pi/om],f)
    y(i)=z($)
  end
endfunction

xmax=2;
x0=.01:.01:xmax; pp=p(x0); 
plot(x0,pp)
//plot2d(x0,pp,rect=[0,0,xmax,xmax])
plot([0,xmax],[0,xmax])

xset("window",1)
t=0:0.01:50/om;
yy=ode(1,0,t,f);
plot(t,yy)


