First question: why is there an extra one at end? Because of the way you handled the "missing" point at the end:
drvY = [drvY , drvY(end)];
with this, what you get is drvY(1999) = drvY(2000). And then what do you do for the second derivative? It's just drvY(1999) - drvY(2000), which is zero. And for the missing point, you are doing same thing again:
secondDrvY = [secondDrvY, secondDrvY(end)]
So you are guaranteed to get two zeros at the end no matter what your data. Handling the boundaries is always hard when doing finite differences. In a general problem, you may get slightly better results if you use gradient() instead of diff(). gradient() uses a second order (i.e. more accurate) central different for the majority of the points, and automatically handles the boundaries (i.e. the returned vector has the same length as the original vector).
However, if all you want to do is the find the zero crossings, then why are you looking at the second derivative to begin with? Just directly look for points that crossed zero:
indices = [find( y(2:end).*y(1:end-1)<0) find(abs(y)<1e-15)]
I think you'll find that works a lot better.
and to answer the other question, here is how I would plot it:
plot(t,y, 'k.', t(indices), zeros(1,length(indices)), 'd', 'MarkerSize', 10, 'MarkerFaceColor', 'r')
