Question: How to correct the code for Delay Differential Equation?

Good day, all.

Please I want to solve the following delay differential equation:

ODE := diff(y(t), t$2) = (2*(1-y(t-1)^2))*(diff(y(t), t))-y(t)

ics := y(0) = 1, (D(y))(0) = 0

using the following codes but there is an error. Please kindly help to modify the codes.

restart;
Digits:=30:

f:=proc(n)
	2*(1-(y[n-1])^2)*delta[n]+y[n]:
end proc:

g:=proc(n)
	-4*y[n-1]*delta[n-1]+2*(1-(y[n-1])^2)*f(n)-delta[n]:
end proc:


e1:=y[n+2] = -y[n]+2*y[n+1]+(1/120)*h^2*(-3*h*g(n+2)+3*g(n)*h+16*f(n+2)+16*f(n)+88*f(n+1)):
e2:=h*delta[n] = -y[n]+y[n+1]-(1/1680)*h^2*(-128*h*g(n+1)-11*h*g(n+2)+59*g(n)*h+40*f(n+2)+520*f(n)+280*f(n+1)):
e3:=h*delta[n+1] = -y[n]+y[n+1]+(1/1680)*h^2*(-152*h*g(n+1)-10*h*g(n+2)+32*g(n)*h+37*f(n+2)+187*f(n)+616*f(n+1)):
e4:=h*delta[n+2] = -y[n]+y[n+1]+(1/1680)*h^2*(128*h*g(n+1)-101*h*g(n+2)+53*g(n)*h+744*f(n+2)+264*f(n)+1512*f(n+1)):

inx:=0:
ind:=0:
iny:=1:
h:=1/2:
n:=1:
omega:=10:
u:=omega*h:
N:=solve(h*p = 10, p):

err := Vector(round(N)):
exy_lst := Vector(round(N)):
numerical_y1:=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+1],delta[n+2]:

step := [seq](eval(x, x=c*h), c=1..N):
printf("%6s%45s%45s\n", 
	"h","Num.y","Num.z");
#eval(<vars>, solve({e||(1..4)},{vars}));


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({e||(1..4)},[par1,par2]), {vars}));

	for i from 1 to 2 do
		printf("%6.5f%45.30f%45.30f\n", 
		h*c,res[i],res[i+2]):
		
		numerical_y1[c] := res[i]:
		
		c:=c+1:
	end do:
	iny:=res[2]:
	ind:=res[4]:
	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[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,blue],symbolsize = 20, legend = ["TFIBF"]);

 Thank you, and best regards.

Please Wait...