How can I change my code so that the c matrix is not overwritten

This is my code modelling advection diffusion to obtain circular diffusion graphs, i get the error Index in position 3 exceeds array bounds (must not exceed 1). I understand I get this because for line [C,h]=contourf(X,Y,Z,300); I am changing the matrix into 2d, and therfore the 3rd dimension is non existent. How can I change the code so that I dont have to overwrite the matrix or that the graphs plot at the same time before the new matrix is assigned.

clf
clear
clc
% Inputs
D=1.04e-6; % Thermal diffusivity, m2/s, eg 0.01
r_in=0;
r_out=0.05;
m=35; %number of divided sections between r_out and r_in
delta_r=(r_out-r_in)/m; %section length, m
nr=m+1; %total no. of radial points
nphi=36;%no. of angle  steps
delta_phi=2*pi/nphi; %angle step
t=40000; % total time, s 
u=delta_r/t; % veloctiy in x direction, m/s
nt=450000; % total no. of time steps
delta_t=t/nt; % timestep, s
Crin=0; % Concentration at r_in, ppm
Crout=1000; % Concentration at r_out, ppm
Cin=0; %initial concentration at r_in<r<r_out, ppm
Cmax=max([Crin,Crout]);
% Solution
for i=1:nr
    r(i)=r_in+(i-1)*delta_r;
end                       
r;
disp(u)
phideg=[0:360/nphi:360+360/nphi];
phi=phideg.*pi/180;
c=u*delta_t/delta_r; % Courant Number
d=D*delta_t/delta_r^2; % diffusion number
d1=D*delta_t/(2*delta_r);
Pe=u*r/D; % Peclet Number
%fprintf('\n')
if d<=0.5
    fprintf('solution stable\nd=%10.7f\n', d)
else
    fprintf('solution unstable\nd=%10.7f\n', d)
end
fprintf('\n') 
fprintf('\n')
if c^2<=2*d
    fprintf('solution stable\n c^2 (%4.2f)<=2*d(%4.2f)\n', c^2, 2*d) 
else
    fprintf('solution unstable\n c^2 (%4.2f)>(%4.2f)\n', c^2, 2*d)
end
% Solution
for i=1:nr
    r(i)=r_in+(i-1)*delta_r;
end                       
% Creating initial boundary conditions 
C=zeros(nr,nphi+2,nt);
for k = 1:nt+1
    for j = 1:nphi+2
        for i=1:nr
            if(i==1)
                C(i,j,k)=Crin;
            elseif(i==nr)
                C(i,j,k)=Crout;
            else
                C(i,j,k)=Cin;
            end
        end
    end
end
C;
% Creating C matrix
for k=1:nt
    for j=1:nphi+2
        for i= 2:nr-1
        C(i,j,k+1)= C(i,j,k) +(c/2)*(C(i+1,j,k)-C(i-1,j,k)) + d*(C(i+1,j,k) - 2*C(i,j,k) + C(i-1,j,k))+(d1/r(i))*(C(i+1,j,k)-C(i-1,j,k));
        end
    end
end
C(:,:,1);
C(:,:,nt+1);
C;
C(:,1,1);
C(:,1,nt+1);
C(:,1,:);
%Initial Temperature profile
C1=C(:,:,1);
subplot(2,2,1)
r=r;
[Phi,R]=meshgrid(phi,r);
Z=C1; %C1.*Phi./Phi;
[X,Y,Z]=pol2cart(Phi,R,Z);
%figure(1)
[C,h]=contourf(X,Y,Z,300);
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel(['r in m'])
%final C
C2=C(:,:,nt+1);      %    <-------------- Where the error appears
subplot(2,2,3)
r=r;
[Phi,R]=meshgrid(phi,r);
Z=C2.*Phi./Phi; 
[X,Y,Z]=pol2cart(Phi,R,Z);
%figure(1)
[C,h]=contourf(X,Y,Z,300);
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel(['r in m'])

[Torsten]

Why not naming the matrix in the call to contourf differently (e.g. M) ?

Why do you need the contour matrix at all ?

Answer 1

[Maneet Kaur Bagga]

Hi,

As per my understanding the error is encountered when you are trying to visualize the results of a 3D matrix "C" as "2D" plots at different time steps. The error is because you're trying to access a slice of the matrix "C" that does not exist in a 2D context after transforming it via "pol2cart". The transformation and contour plotting overwrite the "C" variable.

To resolve this, you should avoid overwriting the variable C with the contour plot handle. In MATLAB, [C,h]=contourf(...) assigns the contour matrix to C, which is also the name of your data matrix. This overwrites your original C matrix when you meant to only create a contour plot.

The possible workaround for this is:

Change the variable names for contour plots: Use a different variable name for the contour plot handles to avoid overwriting the "C" matrix. For example, use [Ch,h] instead of [C,h] when calling contourf.

Please refer to the code snippet below for reference:

% Initial Temperature profile
C1 = C(:,:,1);
subplot(2,2,1)
[Phi,R] = meshgrid(phi,r);
Z = C1; % No need to modify Z here
[X,Y,Z] = pol2cart(Phi,R,Z);
[Ch,h] = contourf(X,Y,Z,300); % Changed C to Ch to avoid overwriting
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Initial Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel('r in m')

% Final C
C2 = C(:,:,nt+1);
subplot(2,2,3)
[Phi,R] = meshgrid(phi,r);
Z = C2; % Directly use C2, assuming you want to plot it as is
[X,Y,Z] = pol2cart(Phi,R,Z);
[Ch,h] = contourf(X,Y,Z,300); % Changed C to Ch to avoid overwriting
set(h,'LineColor','none')
grid on
axis equal
colorbar
title({'Final Concentration Profile - 1D - Circular coordinates';'C in ppm'})
xlabel('r in m')

Hope this helps!