Need help with this MATLAB HW

%% Steady-state diffusion equation
%% PDE and boundary conditions
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
clc; clear;
 
%% add paths
addpath('FVTool');
FVToolStartUp
 
%% Define the domain and create a mesh structure
L = 1;  % domain length
Nx = 20; % number of cells
dx = L/Nx;
m = createMesh2D(Nx,Nx, L,L);
 
x = m.cellcenters.x;
y = m.cellcenters.y;
 
%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=10; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
BC.top.a(:) = 1; BC.top.b(:)=0; BC.top.c(:)=bc(x,1); % top boundary
BC.bottom.a(:) = 1; BC.bottom.b(:)=0; BC.bottom.c(:)=bc(x,0); % bottom boundary
 
 
%% define the transfer coeffs
D_x = 1;
D_y = 1;
D = createFaceVariable(m, [D_x, D_y]);
Mdiff = diffusionTerm(D);
 
%% define source term
 
[X, Y] = ndgrid(x, y);
S = @(X,Y)(6*X.*Y.*(1-Y)-2*X.^3)
s1 = constantSourceTerm(createCellVariable(m,S(X,Y)));
%s1 = 0*s1;
 
[Mbc, RHSbc] = boundaryCondition(BC);
 
M = Mdiff+Mbc;
RHS = -s1+RHSbc;
c = solvePDE(m,M, RHS);
 
figure(1);visualizeCells(c);caxis([0,10]);drawnow;
 
 
% view 3D plot
% add ghost points to x and y
x = [-0.5*dx; x; L+0.5*dx];
y = [-0.5*dx; y; L+0.5*dx];
[X, Y] = ndgrid(x, y);
 
% fix corner points as ghost point
u = c.value;
u(1,1) = u(1,2);
u(1,end) = u(1,end-1);
u(end,1) = u(end,2);
u(end,end) = u(end,end-1);
 
figure(2); surf(X, Y, u);
%axis([0 L 0 L 0 10]);
 
 
function g = bc(x,y)
 
g = -sin(10*pi*x);
 
end

[Walter Roberson]

What difficulty are you asking for assistance with?