You're on the right track trying to use stepinfo. The stepinfo function returns a structure containing all the common step response metrics. You can directly plot the these results. Here is some example code to get you started:
% Defining an arbitrary system
s = tf('s');
G = 1/(s^2 + 5*s + 2);
% Computing the system step response
[y,t] = step(G);
stepResults = stepinfo(y,t);
% Extract the settling time from the stepinfo structure
tSettle= stepResults.SettlingTime;
% Other step response results of interest can be found by looking in
% the stepResults structure
% Find the first index where time exceeds the settling time
% To further improve this, you could interpolate between the points.
indexSettle = find(t >= tSettle,1,'first');
xSettle = t(indexSettle);
ySettle = y(indexSettle);
% Plot step response
plot(t,y);
grid on;
hold on;
% Plot settling time point
plot(xSettle,ySettle,'ro');
% Plot chart lines to settling point
plot([0 xSettle],[ySettle ySettle],'g-.');
plot([xSettle xSettle],[0 ySettle],'g-.');
hold off;
You could easily repeat this same process for each step response metric of interest. Hope this helps!
