% m58_allencahn.m - Allen-Cahn eq. u_t = eps*u_xx+u-u^3, u(-1)=-1, u(1)=1
%         t: forward Euler   x: Chebyshev spectral

% Calling polyval(polyfit(..)) as in this code is exponentially
% unstable, and only works because N is as small as 24.
% One should really use the barycentric formula.

% Differentiation matrix and initial data:
  warning off
  N = 30;
  [D,x] = cheb(N);
  D2 = D^2;
  k = .0025;
  v = .53*x + .47*sin(-1.5*pi*x); vold = v;
  hold off, shg
  xx = linspace(-1,1);
  vv = polyval(polyfit(x,v,N),xx);       % unstable!  should use barycentric
  plt = plot(xx,vv,'.-','linewidth',4);  % polynomial interpolant
  axis([-1 1 -2 2]), grid on
  eps = input('eps? (try .02, .01, .008, .007, .0067, ...) ');

% Solve PDE and stop when in steady state:
  n = 0;
  while 1
    n = n+1;
    v = v + k*(eps*D2*v + v - v.^3);
    v([1 end]) = [1 -1];
    if mod(n,50)==0
       vv = polyval(polyfit(x,v,N),xx); 
       set(plt,'ydata',vv), drawnow    
       title(sprintf('t = %6.2f', n*k), 'fontsize',12)
       if all(diff(v)<=.01) & norm(v-vold)<.001, break, end
       vold = v;
    end
  end