1. To get continuous fit lines, Instead of
plot(x1,y1_fit,'r-')
do something like
x = min(x1):max(x1);
y = polyval(p,x);
plot(x,y,'r--')
2. Have you tried fitting all the data (from 1 to 24) as a single set, using a lower order polynomial. Looks almost quadratic or cubic taken as a whole.
3. You can use the set command to set the MarkerIndices of the x-axis.