clear all
close all
%------------------------------
% Solution of the 1d FKPP equation
% u_t=D*u_xx+u-u^2
% D-diffusion coefficient;
% Initial condition: u(0,x)=f(x);
% x in (-L,L) with
% no-flux boundary conditions: u_x(-L,x)=0, u_x(L,x)=0.
%-----------------------------
% The scheme: the BTCS implicit method 
% u^n+1=A^-1*b
% a=D*dt/dx^2;
%-----------------------------
%-----------------------------
L=50; % length of the string
D=1;
M=500; % amount of points on the (0,L)
x=linspace(-L, L, M+1);
dx=x(2)-x(1);% space grid size
dt=0.05; % time grid size
a=(D*dt)/dx^2;
T=800; % amount of time steps
tt=100;
ys=1.1;
rad=0.01;
nn=1;
a1=-a;
a2=1+2*a;
a3=-a;
%--------------------------------------
u=zeros(T+1,M+1);


f=inline('exp(-5*y.^2)/20','y'); % an initial wave u(x,)

u(1,:)=f(x);

plot(x,u(1,:),'k','LineWidth',2)
hold on
       plot(x,1,'b.','LineWidth',2);
       hold on
       plot(x,0,'r.','LineWidth',2);
axis([-L-rad L+rad 0 ys])
       title(['t=' num2str(0)],'Color','b')
       xlabel('x','FontName','Times','FontSize',10)
       ylabel('u(x,t)','FontName','Times','FontSize',10)
       %saveas(gcf,['FKPP_' num2str(0) '.png']);
       pause(nn)
       hold off
%-----------------------------
%-----------------------------
% Calculate the matrix A, size(A)=M-1xM-1
ah=a2*ones(M+1,1);
ao=a3*ones(M,1);
au=a1*ones(M,1);
A=diag(ah)+diag(au,-1)+diag(ao,1);
A(1,2)=-2*a;
A(M+1,M)=-2*a;
%-----------------------------
%Calculate all other time rows
for t=1:T
    pru=u(t,:)'+dt*(u(t,:)'-u(t,:)'.^2);
    u(t+1,:)=inv(A)*pru;
    if mod(t,tt)==0
       plot(x,u(t+1,:),'k','LineWidth',2);
       hold on
       plot(x,1,'b.','LineWidth',2);
       hold on
       plot(x,0,'r.','LineWidth',2);
       axis([-L-rad L+rad 0 ys])
       title(['t=' num2str(t)],'Color','b')
       xlabel('x','FontName','Times','FontSize',10)
       ylabel('u(x,t)','FontName','Times','FontSize',10)
       hold on 
       %saveas(gcf,['FKPP_' num2str(t) '.png']);
       pause(nn)
       hold off
    end;
end;

%tz=0:T;
% figure(2)
% [C,h]=contourf(x,tz,u,100);
% axis square
% set(h,'LineStyle','none')
% xlabel('x')
% ylabel('t')
% colorbar
% grid off
% 
% figure(3)
% h=surf(x,tz,u)
% shading interp
% light('Position',[-2,2,20])
% lighting phong
% material([0.4,0.6,0.5,30])
% set(h,'EdgeColor','none','FaceColor',[0.2,0.2,0.9],'BackFaceLighting','lit')
% alpha(0.8)
% xlabel('x')
% ylabel('t')
% zlabel('u(x,t)')
% axis([-L-rad L+rad 0 T 0 1])
