clear;
clc;

xl=0;
xr=3*pi; 
tfinal=30;
dx=0.1*pi;
dt=0.01;

nx=(xr-xl)/dx;
nt=tfinal/dt+1;

d1=0.02;
d2=1;
a=1;
b=2;
fsize=15;

s=d1*dt/(dx*dx);
r=d2*dt/(dx*dx);

xmesh=linspace(xl,xr,nx);
tmesh=linspace(0,tfinal,nt);

U=zeros(nx,nt);
V=zeros(nx,nt);
U(:,1)=0.5+0.1*sin(xmesh');
V(:,1)=1.366+0.1*sin(xmesh');

A=zeros(nx,nx);
A(1,1)=1+2*s;
A(1,2)=-2*s;
A(1,nx)=0;
A(nx,1)=0;
A(nx,nx-1)=-2*s;
A(nx,nx)=1+2*s;
for i=2:nx-1
    A(i,i)=1+2*s;
    A(i,i-1)=-s;
    A(i,i+1)=-s;
end

AA=inv(A);

B=zeros(nx,nx);
B(1,1)=1+2*r;
B(1,2)=-2*r;
B(1,nx)=0;
B(nx,1)=0;
B(nx,nx-1)=-2*r;
B(nx,nx)=1+2*r;
for i=2:nx-1
    B(i,i)=1+2*r;
    B(i,i-1)=-r;
    B(i,i+1)=-r;
end
 
BB=inv(B);

for i=1:nt-1   %%for time t
    for j=1:nx %%for space x
        F(j)=U(j,i)+dt*(a-(b+1)*U(j,i)+U(j,i)^2*V(j,i));
        G(j)=V(j,i)+dt*(b*U(j,i)-U(j,i)^2*V(j,i));
    end
    U(:,i+1)=AA*(F');
    V(:,i+1)=BB*(G');
end

surf(tmesh,xmesh,U,'EdgeColor','none');
colormap Jet
colorbar
ylabel('Space x','fontsize',fsize)
xlabel('Time t','fontsize',fsize)
zlabel('u concentration','fontsize',fsize)
view([-1 -1 1])
%view([0 0 1])
figure;
surf(tmesh,xmesh,V,'EdgeColor','none');
colormap Jet
colorbar
ylabel('Space x','fontsize',fsize)
xlabel('Time t','fontsize',fsize)
zlabel('v concentration','fontsize',fsize)
view([-1 -1 1])
%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])
figure;
surf(tmesh,xmesh,V,'EdgeColor','none');
colormap Jet
colorbar
ylabel('Space x','fontsize',fsize)
xlabel('Time t','fontsize',fsize)
view([0 0 1])