PDE convection-diffusion equation using the fully implicit method.

I have wrote thd code below. While runing it I got the following errores that I need help to fix them:

Operator '+' is not supported for operands of type 'struct'.

Error in HW7>createMesh3D (line 394)

G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);

Error in HW7>createMesh3D (line 408)

MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);

Error in HW7 (line 350)

m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh

394 G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);

The code:

%CP7
%Written by Viridiana Salazar
clear variables
close all
%% 1. Space discretization
Lx = 1.0; dx =0.125; N=Lx/dx+1; x=0:dx:Lx;
Ly = 1.0; dy =0.025; y=0:dy:Ly;
Lz = 1.0; dz =0.025; z=0:dz:Lz;
%% 2. Time discretization
tf =0.5; dt =0.001; M=tf/dt+1; t= 0:dt:tf;
%% Constants
Mu=0.5; r=Mu*dt/(dx)^2; q=-0.5*dt/dx;
%% Analytical Solution (We're gonna take from here the initial and
% boundary conditions)
UA = zeros(N,N,N,M);
for n=1:M %time
   for i=1:N %x
      for j=1:N %y
          for k=1:N %z
              UA(i,j,k,n)=(1+exp(x(i)/(2*Mu)+y(j)/(2*Mu)+z(k)/(2*Mu)-3*t(n)/(4*Mu)))^(-1);
          end
      end
   end
end
%% Boundary Conditions
U = zeros(N,N,N,M);
% For x=1 and x=N
U(1,:,:,:)=UA(1,:,:,:); U(N,:,:,:)=UA(N,:,:,:);
%For y=1 and y=N
U(:,1,:,:)=UA(:,1,:,:); U(:,N,:,:)=UA(:,N,:,:);
%For z=1 and z=N
U(:,:,1,:)=UA(:,:,1,:); U(:,:,N,:)=UA(:,:,N,:);
%% Initial conditions
% For t=0
U(:,:,:,1)=UA(:,:,:,1);
%%
for n=1:M-1 %Time loop
% Initial guest at level n+1
for i=1:N %x
   for j=1:N %y
       for k=1:N %z
           if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N %At the boundaries 1 and N, U is known
             U(i,j,k,n+1)=U(i,j,k,n);
           end
       end
   end
end
iteration=0;
EPS=1;
while EPS>0.0000000001 && iteration<20
%Jacobian Matrix
J=zeros(N^3,N^3);
% Vector f
F=zeros(N^3,1);
i=0;
for I=1:N^3
    K=I;
% Fix the index i,j,k
    if mod(I,N*N)==1
       i=i+1;
       j=0;
    end
    if mod(I,N)==1
    k=1;
    j=j+1;
    else
    k=k+1;
    end
    if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
       J(K,I-1)=q+r; %DF/U(i,j,k-1)
       J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
       J(K,I+1)=r; %DF/U(i,j,k+1)
       J(K,I+N)=r; %DF/U(i,j+1,k)
       J(K,I-N)=q+r; %DF/U(i,j-1,k)
       J(K,I+N*N)=r; %DF/U(i+1,j,k)
       J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
       F(K,1)=(q*(-3*U(i,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j,k-1,n+1)+U(i,j,k-1,n+1))+(-(6*r+1)*U(i,j,k,n+1)+r*(U(i+1,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j+1,k,n+1)+U(i,j-1,k,n+1)+U(i,j,k+1,n+1)+U(i,j,k-1,n+1)))+U(i,j,k,n));
    end
%% derivatives for i=N and/or j=N and/or k=N
if i==N && j==N && k==N
   J(K,I-1)=q+r; %DF/U(i,j,k-1)
   J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
   J(K,I-N)=q+r; %DF/U(i,j-1,k)
   J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i==N && j~=N && k~=N
if i==N && j~=N && k~=N && j~=1 && k~=1
   J(K,I-1)=q+r; %DF/U(i,j,k-1)
   J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
   J(K,I+1)=r; %DF/U(i,j,k+1)
   J(K,I+N)=r; %DF/U(i,j+1,k)
   J(K,I-N)=q+r; %DF/U(i,j-1,k)
   J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j==1 && k~=1
   J(K,I-1)=q+r; %DF/U(i,j,k-1)
   J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
   J(K,I+1)=r; %DF/U(i,j,k+1)
   J(K,I+N)=r; %DF/U(i,j+1,k)
   J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j~=1 && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j==N && k~=N
if i~=N && j==N && k~=N && i~=1 && k~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k~=N && i==1 && k~=1
    J(K,I-1)=q*+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j==N && k~=N && i~=1 && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j~=N && k==N
if i~=N && j~=N && k==N && i~=1 && j~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j~=N && k==N && i==1 && j~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j~=N && k==N && i~=1 && j==1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j==N && k~=N
if i==N && j==N && k~=N && k~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j==N && k~=N && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j~=N && k==N
if i==N && j~=N && k==N && j~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q*+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k==N && j==1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i~=N && j==N && k==N
if i~=N && j==N && k==N && i~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k==N && i==1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end

%% derivatives for i=1 and/or j=1 and/or k=1
if i==1 && j==1 && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i==1 && j~=1 && k~=1 && j~=N && k~=N
if i==1 && j~=1 && k~=1 && j~=N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j==N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j~=N && k==N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i~=1 && j==1 && k~=1
if i~=1 && j==1 && k~=1 && i~=N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i==N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i~=N && k==N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=1 && j~=1 && k==1
if i~=1 && j~=1 && k==1 && i~=N && j~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i==N && j~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i~=N && j==N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)  
end
%i==1 && j==1 && k~=1
if i==1 && j==1 && k~=1 && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j==1 && k~=1 && k==N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i==1 && j~=1 && k==1 && j~=N
if i==1 && j~=1 && k==1 && j~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k==1 && j==N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i~=1 && j==1 && k==
if i~=1 && j==1 && k==1 && i~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k==1 && i==N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
end
% Solving the system of equations J*dU=-F, we can use LU decomposition to save
% computational resources in robust system
dU=(J\(-F));
EPS=max(abs(dU(:,1)));
a=0;
for i=1:N %x
    for j=1:N %y
        for k=1:N %z
            a=a+1;
            if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
                U(i,j,k,n+1)=U(i,j,k,n+1)+dU(a);
            end
        end
    end
end
iteration=iteration+1;
end
end


%% Plots and subrotines for 3D plotting
m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh
Numerical=CellVariable(m, U(:,:,:,M));
Analytical=CellVariable(m, UA(:,:,:,M));
Error=CellVariable(m,abs(U(:,:,:,M)-UA(:,:,:,M)));
figure
visualizeCells3D(Numerical);
title('Numerical solution at t=0.25')
figure
visualizeCells3D(Analytical);
title('Analytical solution at t=0.25')
figure
visualizeCells3D(Error);
title('Absolute Error at t=0.25')
function visualizeCells3D(phi)
% Modified from 2012-2016 Ali Akbar Eftekhari
phi.value = phi.value(2:end-1,2:end-1,2:end-1);
[X,Y,Z]=meshgrid(phi.domain.cellcenters.y, phi.domain.cellcenters.x, ...
    phi.domain.cellcenters.z);
phi.value(1)=phi.value(1)+eps; % to avoid an strange error for assigning color limits
Sx = [phi.domain.cellcenters.x(1) phi.domain.cellcenters.x(end)];
Sy = [phi.domain.cellcenters.y(1) phi.domain.cellcenters.y(end)];
Sz = [phi.domain.cellcenters.z(1) phi.domain.cellcenters.z(end)];
slice(X,Y,Z, phi.value, Sy, Sx, Sz);
xlabel('[y vlaues]'); % this is correct [matrix not rotated]
ylabel('[x vlaues]'); % this is correct [matrix not rotated]
zlabel('[z vlaues]');
axis equal tight
colorbar
end
function MS = createMesh3D(varargin)
% Modified from Ali Akbar Eftekhari
Nx=varargin{1};
Ny=varargin{2};
Nz=varargin{3};
Width=varargin{4};
Height=varargin{5};
Depth=varargin{6};
% cell size is dx
dx =0.125;
dy = 0.125;
dz = 0.125;
G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
CellSize.x= dx*ones(Nx+2,1);
CellSize.y= dy*ones(Ny+2,1);
CellSize.z= dz*ones(Nz+2,1);
CellLocation.x= [1:Nx]'*dx;
CellLocation.y= [1:Ny]'*dy;
CellLocation.z= [1:Nz]'*dz;
FaceLocation.x= [0:Nx]'*dx;
FaceLocation.y= [0:Ny]'*dy;
FaceLocation.z= [0:Nz]'*dz;
c=G([1,end], [1,end], [1, end]);
e1=G([1, end], [1, end], 2:Nz+1);
e2=G([1, end], 2:Ny+1, [1, end]);
e3=G(2:Nx+1, [1, end], [1, end]);
MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);
end

% Used for 3D plotting
classdef HW7
    % From 2012-2016 Ali Akbar Eftekhari
    properties
        dimension
        dims
        cellsize
        cellcenters
        facecenters
        corners
        edges 
    end
    methods
        function meshVar = Mesh(dimension, dims, cellsize, ...
                cellcenters, facecenters, corners, edges)
            if nargin>0
                meshVar.dimension = dimension;
                meshVar.dims = dims;
                meshVar.cellsize = cellsize;
                meshVar.cellcenters = cellcenters;
                meshVar.facecenters = facecenters;
                meshVar.corners= corners;
                meshVar.edges= edges;
            end
        end
    end
end

[Torsten]

Seems there is still some error in your code.

You'd better focus on the numerical solution method rather than the graphical representation of the results.

%CP7
%Written by Viridiana Salazar
%clear variables
%close all
%% 1. Space discretization
Lx = 1.0; dx =0.125; N=Lx/dx+1; x=0:dx:Lx;
Ly = 1.0; dy =0.025; y=0:dy:Ly;
Lz = 1.0; dz =0.025; z=0:dz:Lz;

%% 2. Time discretization
tf =0.5; dt =0.001; M=tf/dt+1; t= 0:dt:tf;
%% Constants
Mu=0.5; r=Mu*dt/(dx)^2; q=-0.5*dt/dx;
%% Analytical Solution (We're gonna take from here the initial and
% boundary conditions)
UA = zeros(N,N,N,M);
for n=1:M %time
    for i=1:N %x
        for j=1:N %y
            for k=1:N %z
                UA(i,j,k,n)=(1+exp(x(i)/(2*Mu)+y(j)/(2*Mu)+z(k)/(2*Mu)-3*t(n)/(4*Mu)))^(-1);
            end
        end
    end
end
%% Boundary Conditions
U = zeros(N,N,N,M);
% For x=1 and x=N
U(1,:,:,:)=UA(1,:,:,:); U(N,:,:,:)=UA(N,:,:,:);
%For y=1 and y=N
U(:,1,:,:)=UA(:,1,:,:); U(:,N,:,:)=UA(:,N,:,:);
%For z=1 and z=N
U(:,:,1,:)=UA(:,:,1,:); U(:,:,N,:)=UA(:,:,N,:);
%% Initial conditions
% For t=0
U(:,:,:,1)=UA(:,:,:,1);
%%
for n=1:M-1 %Time loop
% Initial guest at level n+1
UAn = reshape(UA(:,:,:,n),[N,N,N]);
Un = reshape(U(:,:,:,n),[N,N,N]);
disp('Error between analytical and numerical solution:')
norm(UAn(:)-Un(:))
for i=1:N %x
   for j=1:N %y
       for k=1:N %z
           if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N %At the boundaries 1 and N, U is known
             U(i,j,k,n+1)=U(i,j,k,n);
           end
       end
   end
end
iteration=0;
EPS=1;
while EPS>0.0000000001 && iteration<20
%Jacobian Matrix
J=zeros(N^3,N^3);
% Vector f
F=zeros(N^3,1);
i=0;
for I=1:N^3
    K=I;
% Fix the index i,j,k
    if mod(I,N*N)==1
       i=i+1;
       j=0;
    end
    if mod(I,N)==1
    k=1;
    j=j+1;
    else
    k=k+1;
    end
    if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
       J(K,I-1)=q+r; %DF/U(i,j,k-1)
       J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
       J(K,I+1)=r; %DF/U(i,j,k+1)
       J(K,I+N)=r; %DF/U(i,j+1,k)
       J(K,I-N)=q+r; %DF/U(i,j-1,k)
       J(K,I+N*N)=r; %DF/U(i+1,j,k)
       J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
       F(K,1)=(q*(-3*U(i,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j,k-1,n+1)+U(i,j,k-1,n+1))+(-(6*r+1)*U(i,j,k,n+1)+r*(U(i+1,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j+1,k,n+1)+U(i,j-1,k,n+1)+U(i,j,k+1,n+1)+U(i,j,k-1,n+1)))+U(i,j,k,n));
    end
%% derivatives for i=N and/or j=N and/or k=N
if i==N && j==N && k==N
   J(K,I-1)=q+r; %DF/U(i,j,k-1)
   J(K,I)=-3*q-(6*r+1);
%DF/U(i,j,k)
   J(K,I-N)=q+r; %DF/U(i,j-1,k)
   J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i==N && j~=N && k~=N
if i==N && j~=N && k~=N && j~=1 && k~=1
   J(K,I-1)=q+r; %DF/U(i,j,k-1)
   J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)
   J(K,I+1)=r; %DF/U(i,j,k+1)
   J(K,I+N)=r; %DF/U(i,j+1,k)
   J(K,I-N)=q+r; %DF/U(i,j-1,k)
   J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j==1 && k~=1
   J(K,I-1)=q+r; %DF/U(i,j,k-1)
   J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
   J(K,I+1)=r; %DF/U(i,j,k+1)
   J(K,I+N)=r; %DF/U(i,j+1,k)
   J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k~=N && j~=1 && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j==N && k~=N
if i~=N && j==N && k~=N && i~=1 && k~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1);
    %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k~=N && i==1 && k~=1
    J(K,I-1)=q*+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j==N && k~=N && i~=1 && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=N && j~=N && k==N
if i~=N && j~=N && k==N && i~=1 && j~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1);
    %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j~=N && k==N && i==1 && j~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i~=N && j~=N && k==N && i~=1 && j==1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j==N && k~=N
if i==N && j==N && k~=N && k~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1);
    %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i==N && j==N && k~=N && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i==N && j~=N && k==N
if i==N && j~=N && k==N && j~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q*+r; %DF/U(i-1,j,k)
end
if i==N && j~=N && k==N && j==1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
% i~=N && j==N && k==N
if i~=N && j==N && k==N && i~=1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=N && j==N && k==N && i==1
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%% derivatives for i=1 and/or j=1 and/or k=1
if i==1 && j==1 && k==1
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i==1 && j~=1 && k~=1 && j~=N && k~=N
if i==1 && j~=1 && k~=1 && j~=N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j==N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k~=1 && j~=N && k==N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
% i~=1 && j==1 && k~=1
if i~=1 && j==1 && k~=1 && i~=N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i==N && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k~=1 && i~=N && k==N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
%i~=1 && j~=1 && k==1
if i~=1 && j~=1 && k==1 && i~=N && j~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i==N && j~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j~=1 && k==1 && i~=N && j==N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k) 
end
%i==1 && j==1 && k~=1
if i==1 && j==1 && k~=1 && k~=N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j==1 && k~=1 && k==N
    J(K,I-1)=q+r; %DF/U(i,j,k-1)
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i==1 && j~=1 && k==1 && j~=N
if i==1 && j~=1 && k==1 && j~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
if i==1 && j~=1 && k==1 && j==N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I-N)=q+r; %DF/U(i,j-1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
end
%i~=1 && j==1 && k==
if i~=1 && j==1 && k==1 && i~=N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I+N*N)=r; %DF/U(i+1,j,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
if i~=1 && j==1 && k==1 && i==N
    J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)
    J(K,I+1)=r; %DF/U(i,j,k+1)
    J(K,I+N)=r; %DF/U(i,j+1,k)
    J(K,I-N*N)=q+r; %DF/U(i-1,j,k)
end
end
% Solving the system of equations J*dU=-F, we can use LU decomposition to save
% computational resources in robust system
dU=(J\(-F));
EPS=max(abs(dU(:,1)));
a=0;
for i=1:N %x
    for j=1:N %y
        for k=1:N %z
            a=a+1;
            if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N
                U(i,j,k,n+1)=U(i,j,k,n+1)+dU(a);
            end
        end
    end
end
iteration=iteration+1;
end
end

%% Plots and subrotines for 3D plotting
%m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh
%Numerical=CellVariable(m, U(:,:,:,M));
%Analytical=CellVariable(m, UA(:,:,:,M));
%Error=CellVariable(m,abs(U(:,:,:,M)-UA(:,:,:,M)));
%figure
%visualizeCells3D(Numerical);
%title('Numerical solution at t=0.25')
%figure
%visualizeCells3D(Analytical);
%title('Analytical solution at t=0.25')
%figure
%visualizeCells3D(Error);
%title('Absolute Error at t=0.25')
end
function visualizeCells3D(phi)
% Modified from 2012-2016 Ali Akbar Eftekhari
phi.value = phi.value(2:end-1,2:end-1,2:end-1);
[X,Y,Z]=meshgrid(phi.domain.cellcenters.y, phi.domain.cellcenters.x, ...
    phi.domain.cellcenters.z);
phi.value(1)=phi.value(1)+eps; % to avoid an strange error for assigning color limits
Sx = [phi.domain.cellcenters.x(1) phi.domain.cellcenters.x(end)];
Sy = [phi.domain.cellcenters.y(1) phi.domain.cellcenters.y(end)];
Sz = [phi.domain.cellcenters.z(1) phi.domain.cellcenters.z(end)];
slice(X,Y,Z, phi.value, Sy, Sx, Sz);
xlabel('[y vlaues]'); % this is correct [matrix not rotated]
ylabel('[x vlaues]'); % this is correct [matrix not rotated]
zlabel('[z vlaues]');
axis equal tight
colorbar
end
function MS = createMesh3D(varargin)
% Modified from Ali Akbar Eftekhari
Nx=varargin{1};
Ny=varargin{2};
Nz=varargin{3};
Width=varargin{4};
Height=varargin{5};
Depth=varargin{6};
% cell size is dx
dx =0.125;
dy = 0.125;
dz = 0.125;
G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);
CellSize.x= dx*ones(Nx+2,1);
CellSize.y= dy*ones(Ny+2,1);
CellSize.z= dz*ones(Nz+2,1);
CellLocation.x= [1:Nx]'*dx;
CellLocation.y= [1:Ny]'*dy;
CellLocation.z= [1:Nz]'*dz;
FaceLocation.x= [0:Nx]'*dx;
FaceLocation.y= [0:Ny]'*dy;
FaceLocation.z= [0:Nz]'*dz;
c=G([1,end], [1,end], [1, end]);
e1=G([1, end], [1, end], 2:Nz+1);
e2=G([1, end], 2:Ny+1, [1, end]);
e3=G(2:Nx+1, [1, end], [1, end]);
MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);
end
% Used for 3D plotting
%classdef HW7
%    % From 2012-2016 Ali Akbar Eftekhari
%    properties
%        dimension
%        dims
%       cellsize
%        cellcenters
%        facecenters
%        corners
%        edges 
%    end
%    methods
%        function meshVar = Mesh(dimension, dims, cellsize, ...
%                cellcenters, facecenters, corners, edges)
%            if nargin>0
%                meshVar.dimension = dimension;
 %               meshVar.dims = dims;
 %               meshVar.cellsize = cellsize;
 %               meshVar.cellcenters = cellcenters;
 %               meshVar.facecenters = facecenters;
 %               meshVar.corners= corners;
 %               meshVar.edges= edges;
 %           end
 %       end
  %  end
%end

[Hana Bachi]

I listed above the errors I got after running the code. what are the other errors so I can fix them?

[Torsten]

what are the other errors so I can fix them?

I don't know. I just ran your code without graphics and found that the error between analytical solution and numerical approximation grows from timestep to timestep.

For your assignment, the graphical routines of the code are not necessary. According to your description, you should only compute the 2-norm of the difference between numerical and analytical solution at t = 0.5.

[Hana Bachi]

While trying to run the code to get the 3D plots, I'm getting this error message:

Unrecognized function or variable 'CellVariable'.

Error in Cp72nd (line 344)

Numerical=CellVariable(m, U(:,:,:,M));

How can this be fixed?

the code I'm using as a separate .m file is the following

classdef MeshStructure

%CELLVARIABLE Summary of this class goes here

% Detailed explanation goes here

% Copyright (c) 2012-2016 Ali Akbar Eftekhari

% See the license file

properties

dimension

dims

cellsize

cellcenters

facecenters

corners

edges

end

methods

function meshVar = MeshStructure(dimension, dims, cellsize, ...

cellcenters, facecenters, corners, edges)

if nargin>0

meshVar.dimension = dimension;

meshVar.dims = dims;

meshVar.cellsize = cellsize;

meshVar.cellcenters = cellcenters;

meshVar.facecenters = facecenters;

meshVar.corners= corners;

meshVar.edges= edges;

end

end

end

end

[Torsten]

You copied an incomplete code.

The necessary functions are here:

https://de.mathworks.com/matlabcentral/fileexchange/46637-a-simple-finite-volume-solver-for-matlab

But all you need for visualization can be done using MATLAB's "slice" command:

https://de.mathworks.com/help/matlab/ref/slice.html

[kmaonme]

thanks bro! More details for all

https://de.mathworks.com/matlabcentral/fileexchange/46637-a-simple-finite-volume-solver-for-matlab/flagle

[Mary Barber]

By discretizing the equation in both time and space, we can construct a system of equations that can be solved iteratively to obtain the solution at each time step. This method is particularly useful for cases where stability is a concern, but it comes with a higher computational cost https://de.mathworks.com/matlabcentral/fileexchange/46637-a-simple-finite-volume-solver-for-matlab/uno online