Solutions

1 - Setting up Finite Differences

The function FDM should look like:

function M = FDM(N)
%The function FDM returns the Finite Difference Matrix for an N-by-N
%uniform mesh of the unit square.
%INPUTS:
%   N : Number of meshpoints in each direction. Square uniform mesh assumed
%OUTPUTS:
%   M : Finite difference matrix corresponding to the mesh as defined.
% Write the identity matrix of size N to a variable
I = eye(N);
% Write the one-dimensional finite-difference matrix
A = 2*I - diag(ones(N-1,1),1) - diag(ones(N-1,1),-1);
% Create the matrix M
M = (N+1)^2*(kron(A,I) + kron(I,A));
end %function

The script WorksheetScript_Task1_3 should look like:

%% Task 1 - Constructing M
% Write function FDM(N)
% Request input N, construct vectors x and y, then meshgrid X,Y.
N = input('Number of interior points to use (N): ');
x = linspace(0,1,N+2);  x = x(2:end-1); %shave off boundary points
y = x;
[X, Y] = meshgrid(x,y); %create meshgrid of points
M = FDM(N); %create finite difference matrix

Number of interior points to use (N):5

2 - Source and Analytic Functions

The functions F and Analytic should look like this:

function [ val ] = F( z )
%Evaluates the source term of the PDE. 
%INPUTS:
%   z : n x 2 matrix where each row contains a coordinate pair.
%OUTPUTS:
%   val : The source term evaluated at each of the input co-ordinate pairs.
% Calculate the source ensuring that element-wise multiplication is
% used
val = -4*pi*sin(2*pi*z(:,1)).*(cos(2*pi*(z(:,2)-0.5).^2).*pi.*(1+(2*z(:,2)-1).^2)+sin(2*pi*(z(:,2)-0.5).^2));
end %function
function [ val ] = Analytic( z )
%Evaluates the analytic solution of the PDE.
%INPUTS:
%   z : n x 2 matrix where each row contains a coordinate pair.
%OUTPUTS:
%   val : evaluation of the analytic solution at each of the input co-ordinate pairs.
val = sin(2*pi*z(:,1)).*cos(((z(:,2)-0.5).^2)*2*pi);
end %function

3 - Solving the System

At this point your WorksheetScript_Task1_3 should look something like this:

%% Task 0 - Preliminaries
% Good practice to delete all variables and clear all figures before we
% start.
clear;
close all;
%% Task 1 - Constructing M
% Write function FDM(N)
% Request N, construct vectors x and y, then meshgrid X,Y.
N = input('Number of interior points to use (N): ');
x = linspace(0,1,N+2);  x = x(2:end-1); %shave off boundary points
y = x;
[X, Y] = meshgrid(x,y); %create meshgrid of points
M = FDM(N); %create finite difference matrix
%% Task 2 - Source and analytic functions
% Write the following functions:
% Source function F(z)
% Analytic solution function
%% Task 3 - Solving the system
% Evaluate the source and analytic solution on all (x,y) pairs
source = F([X(:),Y(:)]);
uExact = Analytic([X(:),Y(:)]);
% Solve the system AU = F using the backslash operator
uApprox = -M\source;
% Plot the approximate solution on a surface plot
figure()
surf(x,y,reshape(uApprox,N,N))
% Print the error
err = norm(uApprox - uExact);
fprintf('The error is %.7f \n',err)

Number of interior points to use (N):5
The error is 0.2567419 

4 - Error Analysis

Your new script should now look like this:

%% Task 4 - Error Analysis
% Good practive to clear all our variables and figures before we start
clear;
close all;
% Create a range of values for N
NRange = 5:5:75;
hRange = 1./(NRange+1);
% Pre-allocate a vector to contain the errors
errVec = zeros(1,length(NRange));
% Loop over the values of N, completing task 3 for each
for i = 1:length(NRange)
    % Extract the value of N, and calculate the meshgrid
    N = NRange(i);
    x = linspace(0,1,N+2);
    x = x(2:end-1);
    y = x;
    [X,Y] = meshgrid(x,y);
    % Use FDM to construct the M matrix
    M = FDM(N);
    % Evaluate source and analytic solution on all (x,y) pairs
    source = F([X(:),Y(:)]);
    uExact = Analytic([X(:),Y(:)]);
    %Solve the system AU = -F using the backslash operator
    uApprox = -M\source;
    %Store the error
    errVec(i) = norm(uApprox - uExact);
end %for
% Log-log plot of h vs errVec
figure;
loglog(hRange,errVec);
xlabel('$h$','interpreter','latex');
ylabel('Norm of error vector','interpreter','latex');
% Use polyfit to estimate the rate of convergence
% E = C*h^alpha so ln(E) = alpha*ln(h) + ln(C)
% Fit a degree 1 polynomial
pCoeffs = polyfit(log(hRange), log(errVec), 1);
alpha = pCoeffs(1);
% Print out the estimated convergence rate
fprintf('Estimated convergence rate is %.5f \n',alpha);

Estimated convergence rate is 0.99964 

The value of $\alpha$ should be close to $1.0$, and your log-log plot should look linear!