function [xx,yy,C]=test_diffusion_jacobi()

% Numerical method to solve 
% (d/dt)C=[D_{xx}+D_{yy}]C+s(x,y),
% subject to periodic boundary conditions in the x-direction,
% and Neuman boundary conditions at y=0 and y=L_y.
  
% *************************************************************************
% Simulation parameters

Lx=2.d0;
Ly=1.d0;

A0=10.d0;

dx=Lx/(Nx-1);
dy=Ly/(Ny-1);

kx0=2*pi/Lx;
ky0=pi/Ly;

dt=1.d-3;
t_final=0.4d0;
n_timesteps=floor(t_final/dt);
        
% *************************************************************************% Initialise source and concentration fields.
% Initiali
s_source=zeros(Nx,Ny);
C=zeros(Nx,Ny);

xx=0*(1:Nx);
yy=0*(1:Ny);

kx=kx0;
ky=3*ky0;
  
for i=1:Nx
    for j=1:Ny
      x_val=(i-1)*dx;
      y_val=(j-1)*dy;
      
      xx(i)=x_val;
      yy(j)=y_val;
      
      s_source(i,j)=A0*cos(kx*x_val)*cos(ky*y_val);
      C(i,j)=cos(kx0*x_val)*cos(ky0*y_val)+cos(2*kx0*x_val)*cos(ky0*y_val)+cos(kx0*x_val)*cos(4*ky0*y_val);
    end
end

% *************************************************************************  Enter into time loop now
% Enter into time loop now

for iteration_time=1:n_timesteps

    RHS= get_RHS(C,dt,dx,dy,Nx,Ny);
    RHS=RHS+dt*f_source;

    for sor_iteration=1:max_iteration
        C_temp=C;
        C=do_sor_C(RHS,dt,dx,dy,Nx,Ny);

        % Implement Dirichlet conditions at y=0,y=L_y.
        C(:,1)=C(:,2);
        C(:,Ny)=C(:,Ny-1);
    end

    err1= get_diff(C,C_temp,Nx,Ny);
    display(strcat(num2str(iteration_time),': Residual is ', num2str(err1)))
end

end

% *************************************************************************  Enter into time loop now
% *************************************************************************  Enter into time loop now


function RHS=get_RHS(C,dt,dx,dy,Nx,Ny)

ax=dt/(dx*dx);
ay=dt/(dy*dy);

Diffusion=zeros(Nx,Ny);

for j=2:Ny-1
    for i=1:Nx

        if(i==1)
            im1=Nx-1;
        else
            im1=i-1;
        end 

        if(i==Nx)
            ip1=2;
        else
            ip1=i+1;
        end 

        Diffusion(i,j)=ax*( C(ip1,j)+C(im1,j)-2.d0*C(i,j))+ay*( C(i,j+1)+C(i,j-1)-2.d0*C(i,j) );
    end 
end 

RHS=C+0.5d0*Diffusion;

end