why are my MATLAB arrays returned as 840x840 when they should be 840x1?

This is the code I am running to represent the propagation of different masses of pollution in a series of lakes. I then worked out the concentration of pollution in each lake and tried to plot, however my concentration arrays are appearing 840x840 when they should be 840x1. Can anyone spot why this is happening.
this is my code:

% SECTION 1 - Load variables
clear
close all
clc
% Replace the path within quotation marks below to the folder where the
% data is stored in your computer
cd '/Users/brook/Project_Support_Files/Data/Flow Data/'
FL = dir('*.csv');
% Read the spreadsheet with dates and transfom it into an array
t = readtable('Date.xlsx');
t = table2array(t);
% surface areas and volumes

A1 = 8.1925e10;                     % Superior
V01 = 12000*1E9;

A2 = 1.1685e11;                     % Mi-Huron
V02 = (3500+4900)*1E9;

A3 = 2.5404e10;                     % Erie
V03 = 480 * 1E9;

A4 = 1.9121e10;                     % Ontario
V04 = 1640 * 1E9;
% load diversions
% Terms multiplied convert data from m3/s to mm over the lakes surface area
D1 = readmatrix('D1.csv')*1000*60*60*24*30/A1;
D2 = readmatrix('D2.csv')*1000*60*60*24*30/A2;
D3 = readmatrix('D3.csv')*1000*60*60*24*30/A3;
% load connecting flows
F1 = readmatrix('Fo1.csv')*1000*60*60*24*30/A1;
F2 = readmatrix('Fo2.csv')*1000*60*60*24*30/A2;
F3 = readmatrix('Fo3.csv')*1000*60*60*24*30/A3;
F4 = readmatrix('Fo4.csv')*1000*60*60*24*30/A4;
% load evaporation
E1 = readmatrix('E1.csv');
E2 = readmatrix('E2.csv');
E3 = readmatrix('E3.csv');
E4 = readmatrix('E4.csv');
% load precipitation
P1 = readmatrix('P1.csv');
P2 = readmatrix('P2.csv');
P3 = readmatrix('P3.csv');
P4 = readmatrix('P4.csv');
% load runoff
R1 = readmatrix('R1.csv');
R2 = readmatrix('R2.csv');
R3 = readmatrix('R3.csv');
R4 = readmatrix('R4.csv');
% compute ODEs 5.1 to 5.4
% Multiplication by A_i and 1e-3 converts to m3 per month
dV1 = (  R1 + P1 + D1 - F1 - E1) * A1                       * 1E-3;
dV2 = ( (R2 + P2 - D2 - F2 - E2) * A2 +       (F1 * A1) )   * 1E-3;
dV3 = ( (R3 + P3 - D3 - F3 - E3) * A3 +       (F2 * A2) )   * 1E-3;
dV4 = ( (R4 + P4      - F4 - E4) * A4 + (D3 + F3) * A3  )   * 1E-3;
% Integration loops. Timestep = 1month
% Notice that V_i (in this case) is simply the cumulative sum of the
% differentials. For concentration solutions a Euler scheme needs to be
% implemented
V1 = zeros(size(dV1));
V2 = zeros(size(dV2));
V3 = zeros(size(dV3));
V4 = zeros(size(dV4));
for in = 1 : length(dV1)
    V1(in) = sum(dV1(1:in));
    V2(in) = sum(dV2(1:in));
    V3(in) = sum(dV3(1:in));
    V4(in) = sum(dV4(1:in));
end
clear in
V1 = V1 + V01;
V2 = V2 + V02;
V3 = V3 + V03;
V4 = V4 + V04;
%% Plot volume results
% Use this template of plots for your results for concentrations
% Notice that the subplot function helps you plot several graphs in a
% single figure
figure()
subplot(2,2,1)
plot(t,V1*1e-9,'k', 'linewidth', 2)
ylabel('V [km^3]')
title('Superior')
grid on
subplot(2,2,2)
plot(t,V2*1e-9,'b', 'linewidth', 2)
title('Michigan-Huron')
grid on
subplot(2,2,3)
plot(t,V3*1e-9,'c', 'linewidth', 2)
xlabel('Year')
title('Erie')
ylabel(' V [km^3]')
grid on
subplot(2,2,4)
plot(t,V4*1e-9,'r', 'linewidth', 2)
xlabel('Year')
title('Ontario')
grid on
% Concentration of run-off pollution
CR1=0;
CR2=0;
CR3=0;
CR4=0;
%Converting R to litre per month from mm over surface area
R1=R1*A1*(10^(-3))*1000;
R2=R2*A2*(10^(-3))*1000;
R3=R3*A3*(10^(-3))*1000;
R4=R4*A4*(10^(-3))*1000;
%Converting F to litre per month from mm over surface area
F1=F1*A1*(10^(-3))*1000;
F2=F2*A2*(10^(-3))*1000;
F3=F3*A3*(10^(-3))*1000;
F4=F4*A4*(10^(-3))*1000;
%Converting V to litres from m3
V1=V1*1000;
V2=V2*1000;
V3=V3*1000;
V4=V4*1000;
% M values in micrograms (for initial mass of pollution 75000 tonnes )
M01=75000*(10^6)*(10^6);
M02=0;
M03=0;
M04=0;
%Time step
h=1;
t=1:h:length(t);
M1=zeros(size(t));
M2=zeros(size(t));
M3=zeros(size(t));
M4=zeros(size(t));
M1(1)=M01;
M2(1)=M02;
M3(1)=M03;
M4(1)=M04;
% Approximate solution for next values of M1,M2, M3 and M4:
for i=1:(length(t)-1)
dM1=CR1*R1-(M1(i)./V1).*F1;
M1(i+1)=M1(i)+dM1(i)*h;
dM2=CR2*R2+(M1(i)./V1).*F1-(M2(i)./V2).*F2;
M2(i+1)=M2(i)+dM2(i)*h;
dM3=CR3*R3+(M2(i)./V2).*F2-(M3(i)./V3).*F3;
M3(i+1)=M3(i)+dM3(i)*h;
dM4=CR4*R4+(M3(i)./V3).*F3-(M4(i)./V4).*F4;
M4(i+1)=M4(i)+dM4(i)*h;
end
C1=(M1./V1);

C2=(M2./V2);

C3=(M3./V3);

C4=(M4./V4);
% M values in micrograms ( for initial mass of pollutant 50,000 tonnes)
M012=50000*(10^6)*(10^6);
M022=0;
M032=0;
M042=0;
M12=zeros(size(t));
M22=zeros(size(t));
M32=zeros(size(t));
M42=zeros(size(t));
M12(1)=M012;
M22(1)=M022;
M32(1)=M032;
M42(1)=M042;
% Approximate solution for next values of M12, M22, M32 and M42
for i=1:(length(t)-1)
dM12=CR1*R1-(M12(i)./V1).*F1;
M12(i+1)=M12(i)+dM12(i)*h;
dM22=CR2*R2+(M12(i)./V1).*F1-(M22(i)./V2).*F2;
M22(i+1)=M22(i)+dM22(i)*h;
dM32=CR3*R3+(M22(i)./V2).*F2-(M32(i)./V3).*F3;
M32(i+1)=M32(i)+dM32(i)*h;
dM42=CR4*R4+(M32(i)./V3).*F3-(M42(i)./V4).*F4;
M42(i+1)=M42(i)+dM42(i)*h;
end
C12=(M12./V1);
C22=(M22./V2);
C32=(M32./V3);
C42=(M42./V4);
% M values in micrograms (for initial mass pollution 25000 tonnes)
M013=25000*(10^6)*(10^6);
M023=0;
M033=0;
M043=0;
M13=zeros(size(t));
M23=zeros(size(t));
M33=zeros(size(t));
M43=zeros(size(t));
M13(1)=M013;
M23(1)=M023;
M33(1)=M033;
M43(1)=M043;
% Approximate solution for next values of M13, M23, M33 and M43
for i=1:(length(t)-1)
dM13=CR1*R1-(M13(i)./V1).*F1;
M13(i+1)=M13(i)+dM13(i)*h;
dM23=CR2*R2+(M13(i)./V1).*F1-(M23(i)./V2).*F2;
dM23=reshape(dM23',[840,1]);
M23(i+1)=M23(i)+dM23(i)*h;
dM33=CR3*R3+(M23(i)./V2).*F2-(M33(i)./V3).*F3;
dM33=reshape(dM33',[840,1]);
M33(i+1)=M33(i)+dM33(i)*h;
dM43=CR4*R4+(M33(i)./V3).*F3-(M43(i)./V4).*F4;
dM43=reshape(dM43',[840,1]);
M43(i+1)=M43(i)+dM43(i)*h;
end
C13=(M13./V1);
C23=(M23./V2);
C33=(M33./V3);
C43=(M43./V4);
t = readtable('Date.xlsx');
t = table2array(t);
figure()
subplot(2,2,1)
plot(t,C1,'k',t,C12,'y',t,C13, 'c','linewidth',2)
ylabel('Concentration [micrograms per litre]')
title('Superior')
grid on
subplot(2,2,2)
plot(t,C2,'k',t,C22,'y',t,C23, 'c','linewidth',2)
title('Michigan-Huron')
grid on
subplot(2,2,3)
plot(t,C3,'k', t,C32,'y', t,C33, 'c','linewidth',2)
xlabel('Year')
title('Erie')
ylabel(' concentration [micrograms per litre]')
grid on
subplot(2,2,4)
plot(t,C4,'k', t,C42,'y',t,C43,'c','linewidth',2)
xlabel('Year')
title('Ontario')
grid on

I attacthed the code, which was meant to calculate concentration arrays which are 840x1 but they are returned as 840x840 arrays

From your comment:

>M1 is 1x840 and V1 is 840x1

When you do C1=(M1./V1);, you are accidentally causing MATLAB to do implicit expansion and replicate the arrays across the singleton dimensions (of size 1) to output a square matrix which is 840x840.

If you did this instead then C1 would be a column vector:

C1=(M1(:)./V1(:));

Because the (:) operator reshapes the input arrays into columns, so this is now an operation between two columns, regardless whether they were columns or rows originally.

You could also do something like this to be more explicit, and fix all downstream operations:

M1 = M1(:); % equivalent to M1 = reshape(M1,[],1);
V1 = V1(:); % equivalent to V1 = reshape(V1,[],1);
C1 = M1./V1;

More careful construction of M1, V1, etc in the first place to make sure they're always consistently columns (or rows) would be a (subjectively) better fix.