function dy=fpu1(t,y,N,myalpha);

 %edges:
 D(N+1)=y(2)            -2*y(1) +myalpha*((y(2)-y(1))^2-         y(1)^2);
 D(N+N)=y(N-1)          -2*y(N) +myalpha*(       y(N)^2-(y(N)-y(N-1))^2);
 D(1)=y(N+1);
 D(N)=y(N+N);

 %bulk:
 ii=2:N-1;
 D(N+ii)=y(ii+1)+y(ii-1)-2*y(ii)+myalpha*((y(ii+1)-y(ii)).^2-(y(ii)-y(ii-1)).^2);
 D(ii)=y(N+ii);

 dy=D';