这是我运行的代码，用来表示一系列湖泊中不同污染物的传播。然后我计算出每个湖泊中污染物的浓度，并试图绘制，但是我的浓度数组显示为840x840，而它们应该是840x1。有人能看出为什么会发生这种情况吗？这是我的代码：

% 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

我附上了代码，该代码旨在计算浓度阵列（840x1），但它们返回为840x840阵列