Question: How to plot the error of rkf45 without the exact solution

When one equation has no exact solution. I use rkf45 to solve it.

How can I plot the error figures. Can I use the order:s1 := seq([tt[p], Xrre[p] - Xrre[p - 1]], p = 2..N):

restart;
with(plots):
with(LinearAlgebra):
Digits:=15:
w := 500:
k1 := 1:
k2 := sqrt(2):
k3 := -1-sqrt(2):
u1(t):= t/(t + 1)*exp(I * k1 * w * t) + cos(t) * exp(I * k2 * w * t) + exp(I * k3 * w * t):
u2(t):= exp(I * k1 * w * t) + exp(I * k2 * w * t) + t * exp(I * k3 * w * t):
daesys0:={diff(y1(t),t) = I * y2(t) * y2(t) + u1(t), diff(y2(t),t) = I * y1(t) * y1(t) + u2(t), y1(0)= 0, y2(0) = 1};
S:={y1=y1r+I*y1i,y2=y2r+I*y2i};
eval(daesys0,S);
s1:=(evalc@Re)~(%);
s2:=(evalc@Im)~(%%);
res:=dsolve(s1 union s2,numeric, method =rkf45,abserr=10^(-10),relerr=10^(-10),range=0..5,maxfun = 0,output=listprocedure);
Xr,Yr,Xi,Yi:=op(subs(res,[y1r(t),y2r(t),y1i(t),y2i(t)]));

for k from 1 to 1000 do
tt[k] := 0.005 * k;
Xrre[k] := eval(Xr(0.005*k));
Yrre[k] := eval(Yr(0.005*k));
end do:

N := 1000:
s1 := seq([tt[p], Xrre[p] - Xrre[p - 1]], p = 2..N):
s2 := seq([tt[p], Yrre[p] - Yrre[p - 1]], p = 2..N):
plot([s1],font=[Times, roman, 14]);plot([s2],font=[Times, roman, 14]);
#Here we plot the real parts of the error by rkf45