% solve Fisher's equation c_t = D c_xx + a*c*(1-c) with D > 0 and a >= 0
% given constants.  Solve on the interval 0 < x < L, with 
% BC c_x(0,t)=0, c_x(L,t)=0
% and specified initial conditions
%
% Uses Forward-Euler scheme
%
% Use discrete approximations with spatial points xj = (j-1)*dx for
% j=1,...,N+2  where (N+1)*dx=L, i.e. (N+1)=L/dx
%
% We want to find approximations to the solution of the PDE
% at the N+2 points xj for j=1,...,N+2;
%
% define parameters
  D = 0.001;
  a = 0.;
  L = 2.;
%  
  N = input(' enter number of spatial points N ');
  tend = input(' enter final time tend ');
  a  = input(' enter parameter a ');  
  lambda = input(' enter ratio lambda ');
  plotfreq = input(' time interval between plots ');
  tplot = min(plotfreq,tend);
  dx = L/(N+1);
  dt = lambda*(dx*dx)/D; 
  b1 = D*dt/(dx*dx);
  b2 = a*dt;
  x    = dx*[0:1:N+1];
  cold = zeros(N+2,1);
  cnew = zeros(N+2,1);
%
% define initial conditions
%
  for j=1:N+2;
      xx = x(j);
      if xx < L/4 
	 cold(j) = 0.0;
      elseif xx < L/2
	 cold(j) = (xx-L/4)/(L/4);
      elseif xx < 3*L/4	 
	 cold(j) = (3*L/4 - xx)/(L/4);
      else
	 cold(j) = 0.0;
      end
  end
  plot(x,cold,'r*')
  hold on;
%
% march forward in time
%
  t = dt;
  tcount = 0;
  while(t < tend);
    cnew(1    ) = cold(1    ) + b1*(         -2*cold(1)+2*cold(2))+ ...
	          b2*cold(1    ).*(1-cold(1    ));	
    cnew(2:N+1) = cold(2:N+1) + b1*(cold(1:N)-2*cold(2:N+1)+cold(3:N+2))+ ...
	          b2*cold(2:N+1).*(1-cold(2:N+1));
    cnew(N+2  ) = cold(N+2  ) + b1*(2*cold(N+1) - 2*cold(N+2))+...
	          b2*cold(N+2  ).*(1-cold(N+2  ));	          
    tcount = tcount+1;
      if(t > tplot)
         tplot = tplot+plotfreq;
         plot(x,cnew,'b-')
      end;
%%%%    if mod(tcount,20)==0
%%%%       plot(x,cnew,'b-')
%%%%    end
    cold = cnew;
    t = t + dt;
  end;
name = strcat(' Fishers Eq: D=',num2str(D),' a=',num2str(a))
title(name)