You could simply calculate
[x,y] = ode45(@myode,tspan,1,options);
However, have you really looked at your ode? You have, in effect,
dy/dx = x + x.^2.*y + y.^3;
and you start from y = 1.
y blows up extremely rapidly and the ode solvers won't converge after x is slightly greater than 0.4, so there's no way you are going to get to x = 2000.