We can use some of the methods of the object returned by fit to get the pieces we need to display it.
x = (1:10).';
y = x.^2 + 2*x - 3 + randn(size(x));
P = fit(x, y, 'poly2')
The formula method gives us the expression for the fit with the coefficient names.
F = formula(P)
The coeffnames method gives us the coefficient names and the coeffvalues method the coefficient values.
N = coeffnames(P);
V = coeffvalues(P);
Now we can use string operations to replace the coefficient names with the coefficient values (converted to strings.)
formulaWithValues = replace(F, string(N).', string(V))
That + followed by a - for the constant term is a little awkward looking. Let's fix that.
formulaWithValues = replace(formulaWithValues, "+ -", "- ")
Now we could use formulaWithValues to add a legend entry for the fitted curve (or with the text function to put it in the axes itself.)
plot(x, y, 'o', DisplayName = 'raw data')
hold on
xx = 1:0.25:10;
plot(xx, P(xx), DisplayName = formulaWithValues)
legend show
This gets a little more complicated if your fit has problem-dependent parameters, but it wouldn't be that difficult to use probnames and probvalues in addition to coeffnames and coeffvalues.