clear; 
clc;

xl=0; %left endpoint
xr=10; %right endpoint
tfinal=10; %time to integrate
dx=0.2; %spatial grid size
dt=0.01; %time step size

nx=(xr-xl)/dx+1; %number of spatial subintervals 
nt=tfinal/dt+1; %number of time steps

diff=0.1; %diffusion coefficient
a=1; %parameter: growth rate

s=diff*dt/(dx*dx); %diffusion rate

xmesh=linspace(xl,xr,nx); %vector in x direction
tmesh=linspace(0,tfinal,nt); %vector in t direction

U=zeros(nx,nt); %matrix of solution 
U(:,1)=0.4+0.1*sin(xmesh'); %initial condition

A=zeros(nx,nx); %diffusion matrix
A(1,1)=1+s; %first row
A(1,2)=-s;
A(nx,nx-1)=-s; %last row
A(nx,nx)=1+s;
for i=2:nx-1 %other rows
    A(i,i)=1+2*s;
    A(i,i-1)=-s;
    A(i,i+1)=-s;
end

AA=inv(A); %take the inverse matrix

for i=1:nt-1   %%for time t
    for j=1:nx %%for space x
        F(j)=U(j,i)+dt*(a*U(j,i)-U(j,i)^2); %logistic nonlinearity
    end
    U(:,i+1)=AA*(F'); %multiply the inverse matrix
end

fsize=15; %font size for the label

surf(tmesh,xmesh,U,'EdgeColor','none'); %graph in 3D
colormap Jet
colorbar
ylabel('Space x','fontsize',fsize)
xlabel('Time t','fontsize',fsize)
zlabel('Population','fontsize',fsize)
view([-1 -1 1]) %this view gives 3D look
%view([0 0 1])
figure;
surf(tmesh,xmesh,U,'EdgeColor','none');
colormap Jet
colorbar
ylabel('Space x','fontsize',fsize)
xlabel('Time t','fontsize',fsize)
view([0 0 1]) %this view gives 2D look