clear

OPTIONS = [];

rhsfunc = MyODEFunc;
[t,y] = ode45(rhsfunc{2},[0 1000],[0 0], OPTIONS,30,10);
x0 = y(end,:)';

disp('INTEGRATED ODE, x0 = ');
disp(x0)

opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxStepsize',10);
opt=contset(opt,'MaxNumpoints',200);

[x1,v1]=init_EP_EP(@MyODEFunc,x0,[30;6],[1]);
[x,v,s,h,f]=cont(@equilibrium,x1,[],opt);

clf
figure(1)

cpl(x,v,s,[(length(x0)+1) 1]);
axis([0 400 -50 50])