function [t_vec,diff]=create_diff()

global Nx
global Ny
global n_lines

% *************************************************************************
% File structure information

Nx=201;
Ny=101;

n_lines=Nx*Ny;

% *************************************************************************
% Grab all files in folder

files=dir('diffusion_*.dat');
klim=numel(files);

% Compute analytic solution ***********************************************

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

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

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

kx=kx0;
ky=3.d0*ky0;

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

for i=1:Nx
    for j=1:Ny
      xx(i)=(i-1)*dx;
      yy(j)=(j-1)*dy;
      
      C_true(i,j)=(-A0/(kx*kx+ky*ky))*cos(kx*xx(i))*cos(ky*yy(j));
      
    end
end

% *************************************************************************
% Loop over all files from kk=1 up to the last file with index klim.

diff=0*(1:klim);
t_vec=(1:klim)*.2;

for kk=1:klim
    
    % Count the file and output it to the screen as a string message.
    display(strcat('drawing ',num2str(kk),'th file'))
    
    % Grab the name of the kkth file from the cell array.
    filename=files(kk).name;
    
    % Load the data from the kkth file into arrays of double-precision
    % numbers.
    
    % Here, the subfunction my_openfile() is very similar to the one we had
    % before: it opens the file of name "filename" and reads sequentially
    % through the lines of the file, gathering the data into vectors.  The
    % vectors are reshaped into two-dimensional arrays at the end of the
    % operation.
    
    [~,~,C]=my_openfile(filename);
    
    diff(kk)=sqrt(  sum(sum(  (C-C_true).^2))/(Nx*Ny)  );
    
end

end

% *************************************************************************
% *************************************************************************

function [X,Y,C]=my_openfile(filename)

% Global variables - I did not need to pass them explilcitly to the
% subfunction.

global Nx
global Ny
global n_lines

fid=fopen(filename);

X=0*(1:n_lines);
Y=0*(1:n_lines);
C=0*(1:n_lines);

for i=1:n_lines
    c1=fgetl(fid);
    vec_temp=sscanf(c1,'%f');
    x_temp=vec_temp(1);
    y_temp=vec_temp(2);
    C_temp=vec_temp(3);
    
    X(i)=x_temp;
    Y(i)=y_temp;
    C(i)=C_temp;
end

X=reshape(X,Nx,Ny);
Y=reshape(Y,Nx,Ny);
C=reshape(C,Nx,Ny);

fclose(fid);
    
end

% *************************************************************************