Improve the accuracy/convergence of dlyap

I was playing a little with solving the discrete Lyapunov equation (A * X * A' - X + Q = 0) myself when I came across this question, which also posted an answer to. In this answer I did notice that my algorithm scales much worse in computation time that dlyap (which should be O(n^3) I believe), but the accuracy of dlyap gets much worse when the dimension of the matrices gets above 15 by 15.

For example if I use:

n = 15;
A = diag(ones(1,n-1),1);
Q = eye(n);
p = poly(0.1 * ones(1,n));
A(end,:) = -p(end:-1:2);
P = dlyap(A, Q);
accuracy = norm(A * P * A' - P + Q);

This gives the result accuracy = 8.3891e+08, but if I use my method I get accuracy = 5.3357e-15, which is much much much better. And it also seems to imply that dlyap did not converge to the solution. So my question is if it is possible to improve the accuracy/convergence of dlyap, or if this would just be the limitation the algorithm that is used?

Similar results can also be obtained when doing this with lyap instead.

PS: I am not running the latest MATLAB version, namely R2017a, but I assume that dlyap hasn't changed since then.