% 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.
% and specified initial conditions
%
% Uses Crank-Nicolson 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
%
% Uses Crank-Nicolson discretization
%
% We want to find approximations to the solution of the PDE
% at the N+2 points xj for j=1,...,N+2;
%
% define parameters
  cross = 1;
  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 ');
  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 matrices needed for CN method
%
  [ML,MR,R] = matrixdefine(N,lambda);
%
% define initial conditions
%
  epsilon = 0.1;
  if init == 0,
      xleft = 0.5;
      for j=ZERO:LAST;
          xx = x(j);
          c_old(j) = 1.0 - .50*(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))... 	     
	                 - .50*(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;
%
% define boundary condition vector
% (assumes boundary conditions do NOT vary in time)
%
   bcvector = zeros(N+1,1);
   bcvector(1  ) = (lambda/2)*(g0 + g0);
%
% march forward in time
%
  t = dt;
  tcount = 0;
  while(t < tend);
%
% use 4th order Runge-Kutta to solve dc/dt = a*c*(1-c) for one timestep
%
    k1    = a*c_old(ONE:LAST,1).*(1-c_old(ONE:LAST,1));
    ctemp = c_old(ONE:LAST,1) + .5*dt*k1;
    k2    = a*ctemp.*(1-ctemp);
    ctemp = c_old(ONE:LAST,1) + .5*dt*k2;
    k3    = a*ctemp.*(1-ctemp);
    ctemp = c_old(ONE:LAST,1) +    dt*k3;
    k4    = a*ctemp.*(1-ctemp);
    ctemp = c_old(ONE:LAST,1) + (1/6)*dt*(k1+2*k2+2*k3+k4);
%
% now compute effect of diffusion
%
    [ctempnew,outputflag]=matsolve(ML,MR,ctemp,R,bcvector);
% set c_new at x=0 to g0    
    c_new = [g0;ctempnew(:,1)];
    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
%
    ncount = ncount+1;
    x_star(ncount) = L;
    for j=LASTM1:-1:ZERO
        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;	
    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-','Linewidth',2)
     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-','Linewidth',2)
     wavespeed = sqrt(4*D*a);
     title(strcat(' x^*(t) for c^* =',num2str(c_star),' Speed=',num2str(speed),...
                     ' sqrt(4aD) =',num2str(wavespeed)),'Fontsize',14)     
  end;