% 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(0,t)=g0, c_x(L,t)=0 -- I will use g0=1.
% IC specified initial conditions
%
% Uses Forward Euler scheme
%
% Use discrete approximations with spatial points xj = j*dx for
% j=0,...,N+1  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=0,...,N+1;
% We have a value at xj for j=0 (x=0), so we solve equations
%    to find the approximation for j=1,...N+1.
%
% define parameters
  cross = 0;
  D = 0.001;
  a = 0.;
  L = 10.;
  g0 = 1;
%  
  N = input(' enter number of spatial points N ');
  tend = input(' enter final time tend ');
  a  = input(' enter parameter a ');  
  lambda = input(' enter ratio lambda ');
  init   = input(' which initial data 0,1,2? ');
  plotfreq = input(' time interval between plots ');
  cross    = input(' type 1 to calculate crossings ');
  if cross == 1
    ncount = 0
    c_star = input(' enter c_star ');
  end;
  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];
  c_old = zeros(N+2,1);
  c_new = zeros(N+2,1);
  ZERO    = 1;
  ONE     = 2;
  TWO     = 3;
  THREE   = 4;
  LASTM1  = N+1;
  LAST    = N+2;
%
% define initial conditions
%
  epsilon = 0.1;
  if init == 0,
      xleft = 0.5;
      for j=ZERO:LAST;
          xx = x(j);
          c_old(j) = 1 - .5*(1.0 + tanh((xx-xleft)/epsilon));
      end
  elseif init == 1,
      xleft = .5;
      xright = 1.5;
      epsilon = 0.1;
      for j=ZERO:LAST;
          xx = x(j);
          c_old(j) = 1.0 - .25*(1.0 + tanh((xx-xleft)/epsilon)) ...
	                 - .25*(1.0 + tanh((xx-xright)/epsilon));
      end    
  elseif init == 2,
      xleft = .5;
      xmid  = 1.0;
      xright = 1.5;
      epsilon = 0.1;    
      for j=ZERO:LAST;
          xx = x(j);
          c_old(j) = 1.0 - .25*(1.0 + tanh((xx-xleft )/epsilon))...
	                 + .25*(1.0 + tanh((xx-xmid  )/epsilon))... 	     
	                 -  .5*(1.0 + tanh((xx-xright)/epsilon)); 
      end        
  end;
%  
  subplot(2,1,1)
  plot(x,c_old,'r*'),axis([0 L 0 1.1])
  hold on;

%
% march forward in time
%
  t = dt;
  tcount = 0;
  while(t < tend);
    c_new(ZERO ) = g0;
    c_new(ONE  ) = c_old(ONE) + b1*(    -2*c_old(ONE)+c_old(TWO))+ ...
	          b2*c_old(ONE).*(1-c_old(ONE)) + b1*g0;	
    c_new(TWO:LASTM1) = c_old(TWO:LASTM1) + b1*(c_old(ONE:LAST-2)-2*c_old(TWO:LASTM1)+c_old(THREE:LAST))+ ...
	          b2*c_old(TWO:LASTM1).*(1-c_old(TWO:LASTM1));
    c_new(LAST     ) = c_old(LAST ) + b1*(2*c_old(LASTM1) - 2*c_old(LAST))+...
	          b2*c_old(LAST).*(1-c_old(LAST));	          
    tcount = tcount+1;
      if(t > tplot)
         tplot = tplot+plotfreq;
         plot(x,c_new,'b-')
      end;
    c_old = c_new;
%
% determine location of point at which c_old = c_star
%
    if (cross == 1)
    ncount = ncount+1;
    x_star(ncount) = L;
    for j=LASTM1:-1:1
        if c_old(j) > c_star
           slope = (c_old(j+1)-c_old(j))/dx;
           x_star(ncount) = (c_star-c_old(j))/slope + (j-1)*dx;
	   break;
        end;	      
    end;	
    end;
    t = t + dt;
  end;
name = strcat(' Fishers Eq: D=',num2str(D),' a=',num2str(a))
title(name,'Fontsize',14)
  if cross == 1,
     subplot(2,1,2)
     times = [1:1:ncount]*dt;
     plot(times,x_star,'k-')
     hold on;
     title(strcat(' x^*(t) for c^* =',num2str(c_star)),'Fontsize',14)
%
% do least-square line fit to last m points of (times,x_star)
%
     m = floor(tcount/10);
     ts = times (ncount-m+1:ncount);
     xs = x_star(ncount-m+1:ncount);
     tavg = sum(ts)/m;
     xavg = sum(xs)/m;
     ttavg = sum(ts.*ts)/m;
     xtavg = sum(xs.*ts)/m;
     slope = (xtavg - xavg*tavg)/(ttavg - tavg*tavg);
     inter = (xavg - slope*tavg);
     speed = slope
     xfit = slope*ts + inter;
     plot(ts,xfit,'g-')
     wavespeed = sqrt(4*D*a);
     title(strcat(' x^*(t) for c^* =',num2str(c_star),' Speed=',num2str(speed),...
                     ' sqrt(4aD) =',num2str(wavespeed)),'Fontsize',14)     
%%%%%%     title(strcat(' x^*(t) for c^* =',num2str(c_star),' Speed=',num2str(speed)),'Fontsize',14)     
  end;