完整曲线：

% 假设阀门瞬间关闭
% 初始数据:

clear
tic
L=3000;          % 管线长度
Hr=70;          % 泵压力
N=10;            % 分段数
NS=N+1;         % 节点数
e=0.001651;     % 壁厚m,0.065''
D=0.00635-2*e;    % 管道内径
K=2.1e+9;       % 流体体积弹性系数
Rho=1000;        % 液体密度kg/m^3
E=2.1e11;       % 弹性模数
G=9.806;        % 重力加速度
f=0.018;        % 摩擦系数
A=pi*D^2/4;     % 管道横截面积
dx=L/N;         % 单位长度
t_max=20;        % 总计算时间
a=sqrt(K/Rho/(1+K*D/E*e)); % 波速
dt=dx/a;        % 单位计算时间
k=1;
t=0;time(k)=t;

B=a/(G*A);
R=f*dx/(2*G*D*A^2);
Q0=0;
H0=0;

for j=1:NS                  
Q(k,j)=0; % 第一行流量
H(k,j)=0; % 第一行压头
end
H(1,1)=Hr;

while t<t_max
   for j=2:N
        CP=H(k,j-1)+B*Q(k,j-1)-R*Q(k,j-1)*abs(Q(k,j-1));    % +dt*Q(k,j-1)/A;    % CP==C-==CR
        CM=H(k,j+1)-B*Q(k,j+1)+R*Q(k,j+1)*abs(Q(k,j+1));    % -dt*Q(k,j+1)/A;    % CM==C+==CL
        H(k+1,j)=(CP+CM)/2;       % 计算压头
        Q(k+1,j)=(H(k+1,j)-CM)/(B);   % 计算流量
         
   end
    
   % 边界值需要单独计算
    CP=H(k,N)+B*Q(k,N)-R*Q(k,N)*abs(Q(k,N));    % +dt*Q(k,N)/A;    % 单独计算CP
    CM=H(k,2)-B*Q(k,2)+R*Q(k,2)*abs(Q(k,2));    % -dt*Q(k,2)/A;    % 单独计算CM
    H(k+1,1)=Hr;                                % 压头第一列
    Q(k+1,1)=(H(k+1,1)-CM)/B;                     % 流量第一列数值，边界条件公式
    Q(k+1,NS)=0;                                % 流量最后一列
    H(k+1,NS)=CP;                               % 压头最后一列数值，边界条件公式
    t=t+dt;
    k=k+1;
    time(k)=t;
    
end
toc

plot(time,H(:,N+1))
% hold on
% plot([0,t_max],[Hr,Hr],'b:') % 选取阀门处压力值绘制曲线图
% hold on
% plot([0,t_max],[0,0],'b:') % 选取阀门处压力值绘制曲线图
title('MOC-阀门处圧力曲线');
xlabel('单位：s');
ylabel('单位：m');