I tried running the following codes but got error
restart;
Digits := 30:
f1 := proc(n)
lamda * y1[n] + (y2[n])^2;
end proc:
f2 := proc(n)
-y2[n];
end proc:
F1 := proc(n)
lamda * f1(n) + 2 * f2(n) * y2[n];
end proc:
F2 := proc(n)
-f2(n);
end proc:
e1 := y1[n+2] = (53/485) * y1[n] - (6/485) * F1(n+2) * h^2 + (6/485) * f1(n) * h +
(128/485) * f1(n+1/2) * h + (384/485) * f1(n+3/2) * h +
(432/485) * f1(n+1) * h + (20/97) * f1(n+2) * h +
(432/485) * y1[n+1] + (512/485) * y1[n+1/2] - (512/485) * y1[n+3/2]:
e2 := h^2 * F1(n+1/2) = -(9/970) * F1(n+2) * h^2 + (503/1940) * f1(n) * h -
(9412/1455) * f1(n+1/2) * h - (682/485) * f1(n+3/2) * h -
(4041/485) * f1(n+1) * h + (83/1164) * f1(n+2) * h +
(3878/1455) * y1[n] + (9054/485) * y1[n+1] -
(14166/485) * y1[n+1/2] + (11458/1455) * y1[n+3/2]:
e3 := h^2 * F1(n+1) = -(28/4365) * F1(n+2) * h^2 + (541/8730) * f1(n) * h +
(25072/13095) * f1(n+1/2) * h - (5968/4365) * f1(n+3/2) * h +
(224/485) * f1(n+1) * h + (269/5238) * f1(n+2) * h +
(7532/13095) * y1[n] - (9476/485) * y1[n+1] +
(4864/485) * y1[n+1/2] + (116992/13095) * y1[n+3/2]:
e4 := h^2 * F1(n+3/2) = -(31/970) * F1(n+2) * h^2 + (671/5820) * f1(n) * h +
(3902/1455) * f1(n+1/2) * h + (12676/1455) * f1(n+3/2) * h +
(5481/485) * f1(n+1) * h + (329/1164) * f1(n+2) * h +
(1502/1455) * y1[n] + (9846/485) * y1[n+1] +
(5526/485) * y1[n+1/2] - (47618/1455) * y1[n+3/2]:
e5 := y2[n+2] = (53/485) * y2[n] - (6/485) * F2(n+2) * h^2 + (6/485) * f2(n) * h +
(128/485) * f2(n+1/2) * h + (384/485) * f2(n+3/2) * h +
(432/485) * f2(n+1) * h + (20/97) * f2(n+2) * h +
(432/485) * y2[n+1] + (512/485) * y2[n+1/2] - (512/485) * y2[n+3/2]:
e6 := h^2 * F2(n+1/2) = -(9/970) * F2(n+2) * h^2 + (503/1940) * f2(n) * h -
(9412/1455) * f2(n+1/2) * h - (682/485) * f2(n+3/2) * h -
(4041/485) * f2(n+1) * h + (83/1164) * f2(n+2) * h +
(3878/1455) * y2[n] + (9054/485) * y2[n+1] -
(14166/485) * y2[n+1/2] + (11458/1455) * y2[n+3/2]:
e7 := h^2 * F2(n+1) = -(28/4365) * F2(n+2) * h^2 + (541/8730) * f2(n) * h +
(25072/13095) * f2(n+1/2) * h - (5968/4365) * f2(n+3/2) * h +
(224/485) * f2(n+1) * h + (269/5238) * f2(n+2) * h +
(7532/13095) * y2[n] - (9476/485) * y2[n+1] +
(4864/485) * y2[n+1/2] + (116992/13095) * y2[n+3/2]:
e8 := h^2 * F2(n+3/2) = -(31/970) * F2(n+2) * h^2 + (671/5820) * f2(n) * h +
(3902/1455) * f2(n+1/2) * h + (12676/1455) * f2(n+3/2) * h +
(5481/485) * f2(n+1) * h + (329/1164) * f2(n+2) * h +
(1502/1455) * y2[n] + (9846/485) * y2[n+1] +
(5526/485) * y2[n+1/2] - (47618/1455) * y2[n+3/2]:
Digits := 30: lamda := 10000:
h := 0.1:
N := solve(h*p = 10/4, p):
n := 0:
exy1 := [seq](eval(-exp(-2*i/2)/(lamda+2)), i = h .. N, h):
exy2 := [seq](eval(exp(-i/2)), i = h .. N, h):
iny1 := 1:
iny2 := 1:
c := 1:
vars := y1[n+1/2], y1[n+1], y1[n+3/2], y1[n+2], y2[n+1/2], y2[n+1], y2[n+3/2], y2[n+2]:
printf("%6s%15s%15s%15s%15s%10s%10s\n",
"h", "numy1", "numy2",
"exy1", "exy2",
"erry1", "erry2");
tolerance := 1e-6:
st := time():
for k from 1 to N do
res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars}), tolerance = tolerance):
for i from 1 to 4 do
printf("%6.5f%15.10f%15.10f%15.8g%15.10g%10.3g%10.3g\n",
h*c, res[i], res[i+4],
exy1[c], exy2[c],
abs(res[i]-exy1[c]), abs(res[i+4]-exy2[c])):
c := c+1:
end do:
iny1 := res[4]:
iny2 := res[8]:
end do:
v := time() - st: