modmid := proc(f, a, ya, h, N)
local t, w, i;
t := a;
w[0] := ya;
w[1] := ya + h*f(t, w[0]);

for i from 1 to N-1 do
	w[i+1] := w[i-1] + 2*h*f(t+i*h, w[i]);
	t := t + h;
end do;
return(w[i])
end proc;