hello friend

please run this code and resend result to me.

thank you.

restart;

with(LinearAlgebra); with(plots);

c[1] := 0; c[2] := 1/2;

a[2, 1] := 1/8;

b[-1] := -1/3; b[1] := 2/3; b[2] := 5/6;

bb[2] := 5/12;

s := 2;

e[1] := Matrix(s+1, 1); e[2] := Matrix(s+1, 1); e[3] := Matrix(s+1, 1); e[4] := Matrix(s+1, 1);

c[0] := 0;

for ii to s+1 do e[1][ii, 1] := 1; e[2][ii, 1] := -1; e[3][ii, 1] := c[ii-1] end do;

e[1][1, 1] := 0; e[2][2, 1] := 0; e[4] := -e[3];

A := Matrix(s+1, s+1);

for i from 3 to s+1 do for j from 2 to i-1 do A[i, j] := a[i-1, j-1] end do; A[i, 1] := A[i, 2] end do;

Id := Matrix(s+1, s+1, shape = identity);

N := Matrix(s+1, s+1);

for i to s+1 do for j to s+1 do N[i, j] := -H*A[i, j]+Id[i, j] end do end do;

Bb := Matrix(s+1, 1);

for i from 3 to s+1 do Bb[i, 1] := bb[i-1] end do;

B := Matrix(s+1, 1);

for i from 2 to s+1 do B[i, 1] := b[i-1] end do; B[1, 1] := b[-1];

Z := Multiply(Transpose(Bb), 1/N);

for i to s+1 do Z[1, i] := H*Z[1, i] end do;

Q := Multiply(Transpose(B), 1/N);

for i to s+1 do Q[1, i] := H*Q[1, i] end do;

Dh := proc (H) options operator, arrow; Matrix([[1+Multiply(Z, e[1]), Multiply(Z, e[2]), 3/2+Multiply(Z, e[3]), -1/2+Multiply(Z, e[4])], [1, 0, 0, 0], [Multiply(Q, e[1]), Multiply(Q, e[2]), 1+Multiply(Q, e[3]), Multiply(Q, e[4])], [0, 0, 1, 0]]) end proc;

k := Matrix(4, 4);

for i to 4 do for j to 4 do if i = j then k[i, j] := ep-Dh(H)[i, j] else k[i, j] := -Dh(H)[i, j] end if end do end do;

sp := unapply(collect(Determinant(k), ep), H, ep);

sol := solve(sp(H, ep), ep, explicit);

sol := simplify(map(allvalues, {sol}));

l := numelems(sol);

inequal({evalc(abs(subs(H = x+I*y, sol[1]))) <= 1, evalc(abs(subs(H = x+I*y, sol[2]))) <= 1, evalc(abs(subs(H = x+I*y, sol[3]))) <= 1, evalc(abs(subs(H = x+I*y, sol[4]))) <= 1}, x = -3 .. 3, y = -3 .. 3, color = "Nautical 1");