N = 100;
h = 2*pi/N;
x = (0:N)*h;
k = 0.5*h^2
t = 0 ;
u = sin(x);
flag=0;

plot(x,u)
drawnow

while(t<2)
  unew = u(2:N) + k*((u(3:(N+1))-2*u(2:N)+u(1:(N-1)))/h^2 + (4+sin(x(2:N))).*(u(3:(N+1))-u(1:(N-1)))/2/h);
  u = [0 unew 0]; 
  t = t+k;
  plot(x,u)
  drawnow
  if flag==0 && max(u)==0
      thet=t;
      flag=1;
  end
end

thet

%N=50  k =  0.0079      thet = 0.9475
%N=100 k =  0.0020      thet = 0.9495
%N=200 k =  4.9348e-04  thet = 0.9495
%N=400 k =  1.2337e-04  thet = 0.9495
