Question: How do I modify the following maple code in a simple form to allow the better values of abs(res[i] - exy) to about 10^-11?

Dear Colleague. 

I am trying to improve the results of abs(res[i] - exy) in the following codes.

restart;
Digits := 30:

# Define the function
f := proc(n)
    -0.5*y[n] + 0.5*sin(x[n] - Pi)
end proc:

# Define equations
e1 := y[n+2] = 2*h*delta[n] + y[n] - h^2*(-2*sin(u)*f(n)*u^2 - 2*sin(u)*f(n+2)*u^2 + 2*sin(2*u)*f(n+1)*u^2 + 2*cos(u)*f(n)*u - 2*cos(u)*f(n+2)*u + 2*cos(2*u)*f(n+1)*u - 2*cos(2*u)*f(n)*u - 2*sin(u)*f(n) + 2*sin(u)*f(n+2) + sin(2*u)*f(n) - sin(2*u)*f(n+2) - 2*f(n+1)*u + 2*f(n+2)*u)/((2*sin(u) - sin(2*u))*u^2):
e2 := y[n+1] = h*delta[n] + y[n] - (1/2)*h^2*(-sin(u)*f(n)*u^2 - sin(u)*f(n+2)*u^2 + sin(2*u)*f(n+1)*u^2 + 2*cos(u)*f(n)*u - 2*cos(u)*f(n+2)*u + 2*cos(2*u)*f(n+1)*u - 2*cos(2*u)*f(n)*u + 4*sin(u)*f(n+1) - 4*sin(u)*f(n) - 2*sin(2*u)*f(n+1) + 2*sin(2*u)*f(n) - 2*f(n+1)*u + 2*f(n+2)*u)/((2*sin(u) - sin(2*u))*u^2):
e3 := h*delta[n+2] = h*delta[n] + h^2*(2*sin(u)*f(n)*u + 2*sin(u)*f(n+2)*u - 2*sin(2*u)*f(n+1)*u - 2*cos(2*u)*f(n+1) + cos(2*u)*f(n) + cos(2*u)*f(n+2) + 2*f(n+1) - f(n) - f(n+2))/(u*(2*sin(u) - sin(2*u))):

with(LinearAlgebra):
epsilon := 10^(-10):
inx := 0:
ind := 1:
iny := 0:
h := 0.01:
n := 0:
omega := 1:
u := omega * h:
tol := 1e-4:
N := solve(h * p = 8 * Pi, p):

err := Vector(round(N)):
exy_lst := Vector(round(N)):

c := 1:
for j from 0 to 2 do
    t[j] := inx + j * h:
end do:

vars := y[n+1], y[n+2], delta[n+2]:

step := [seq(eval(x, x = c * h), c = 1 .. N)]:
printf("%6s%15s%15s%16s%15s%15s%15s\n", "h", "Num.y", "Num.z", "Ex.y", "Ex.z", "Error y", "Error z");

st := time():
for k from 1 to N / 2 do
    par1 := x[0] = t[0], x[1] = t[1], x[2] = t[2]:
    par2 := y[n] = iny, delta[n] = ind:    
    
    res := eval(<vars>, fsolve(eval({e1, e2, e3}, [par1, par2]), {vars}));

    for i from 1 to 2 do
        exy := eval(sin(c * h)):
        exz := eval(cos(c * h)):
        printf("%6.5f%17.9f%15.9f%15.9f%15.9f%13.5g%15.5g\n", h * c, res[i], res[i+1], exy, exz, abs(res[i] - exy), abs(res[i+1] - exz));
        
        err[c] := abs(evalf(res[i] - exy));
        if Norm(err) <= tol then 
            h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
        else 
            break
        end if;
        exy_lst[c] := exy;
        numerical_y1[c] := res[i];
        c := c + 1;
    end do;
    iny := res[2];
    ind := res[3];
    inx := t[2];
    for j from 0 to 2 do
        t[j] := inx + j * h;
    end do;
end do:
v := time() - st;
v / 4;
printf("Maximum error is %.13g\n", max(err));
NFE = evalf((N / 4 * 3) + 1);

# Get array of numerical and exact solutions for y1
numerical_array_y1 := [seq(numerical_y1[i], i = 1 .. N)]:
exact_array_y1 := [seq(exy_lst[i], i = 1 .. N)]:

# Get array of time steps
time_t := [seq(step[i], i = 1 .. N)]:

# Display graphs for y1
with(plots):
numerical_plot_y1 := plot(time_t, numerical_array_y1, style = point, symbol = asterisk, color = blue, symbolsize = 20, legend = ["TFIBF"]);
exact_plot_y1 := plot(time_t, exact_array_y1, style = point, symbol = box, color = red, symbolsize = 20, legend = ["EXACT"]);

display({numerical_plot_y1, exact_plot_y1});
Error_plot_y1 := plot(time_t, err, style = line, symbol = box, tickmarks = [piticks, decimalticks], color = navy, labels = [`h=Pi/8`, typeset(`Absolute Errors`)]);

I am suspecting that I didnt update the new h properly (I may be wrong, though). Please kindly help modify the code to allow the values of abs(res[i] - exy) to about 10^(-11). Thank you and best regards.