function F=funcHeNeu(t,y)
%
%
global m eta1 eta0
mm=m+1;
mmm=m+2;
e=ones(mmm,1);
GHeNeu=gHeNeu(t);
%
A=mm^2*spdiags([e -2*e e],-1:1,mmm,mmm);
A(1,1)=A(1,1)-mm*2*eta0;
A(1,2)=A(1,2)+mm^2;
A(mmm,mmm)=A(mmm,mmm)-mm*2*eta1;
A(mmm,mmm-1)=A(mmm,mmm-1)+mm^2;
%
F=A*y+2*mm*GHeNeu;