L=20;
dx=0.2;
N=int(L/dx)
dt=0.1;
Tfin=600;
M=int(Tfin/dt)
v=0.01;
a=0.25;
eps=0.01;
bet=0.2;
deff('y=f(u)','y=u*(1-u)*(u-a)');
gam=1
I0=0
h0=4;
T0=4;
//
u=zeros(1,N+1)
unew=u;
w=zeros(1,N+1)
wnew=u;
u(1)=h0
//
for i=1:M
  if (i*dt<T0) then
      unew(1)=h0
  else
       unew(1)=0
  end       
  for j=2:N
    unew(j)=u(j)+dt*(v/(dx^2)*(u(j-1)-2*u(j)+u(j+1))+f(u(j))-w(j)+I0) 
    wnew(j)=w(j)+dt*eps*(bet*u(j)-gam*w(j))
  end
  u=unew;w=wnew
  drawlater
  clf()
  plot2d(0:dx:L,u,rect=[0,-1,20,5])
  drawnow
end
clf()
plot2d(0:dx:L,u,rect=[0,-1,20,2])
