一室模型4

药物在体内分布与排除的一室模型与分析
摘要本文讨论了药物动力学中的血药浓度随时间变化规律及确定给药方案的问题。

针对问题一，在只有中心室的条件下，运用微分的数学思想，建立了一次给药时血药浓度关于时间的微分方程模型。

按照快速静脉注射、恒速静脉滴注和口服或肌肉注射3种给药方式，不同的初值对应微分方程不同的解，分别得到3种给药方式下血药浓度随时间变化的表达式。

根据某种药物血药浓度随时间变化的数据，利用最小二乘法进行曲线拟合，解出此药物动力学参数。

运用Matlab 软件分别画出3种情况下的血药浓度曲线图，反映了血药浓度随时间递减或波动递减的趋势。

针对问题二和问题三，在快速静脉注射、恒速静脉滴注和口服或肌肉注射3种多次重复给药方式下，采用问题一中一次给药方式下的血药浓度表达式，运用迭代的方法计算出多次给药方式下血药浓度表达式，利用Matlab软件画出血药浓度曲线的图形。

通过控制变量法，依据相关资料分别赋予时间间隔和给药剂量不同的值，运用Matlab软件画出多组血药浓度随时间变化的曲线图,进行分析、比较。

联系具体实际情况，即可对给药时间间隔与给药剂量的确定提供指导。

本文中建立的为线性微分方程模型，可推广至非线性模型，在现实生活中可运用于农药残留物和血液中酒精浓度的测定。

关键词微分方程模型；曲线拟合；迭代；控制变量法
一、问题重述
随着社会的发展，人们越来越关注医疗保健，这极大地推动了临床医学的发展，而在这一过程中药物疗效的测定显得尤为重要。

根据药物动力学知识知，血药浓度即单位体积内药物的含量的大小会直接影响药物的疗效，浓度太低达不到预期的效果，浓度太高会导致药物中毒、副作用太强和药物浪费。

而药物在机体内的一系列化学反应会影响血药浓度随时间、空间的变化。

药物动力学在研究血药浓度随时间变化的过程中常建立房室模型，房室是指有机体的一部分，机体会因为其研究过程中条件和对象的不同分为多个房室，同一房室内的药物分布均匀，即血药浓度相同。

讨论按固定时间间隔，多次重复按固定剂量给药方式，为了维持药品的疗效和保证机体的安全，要求血药浓度在一个合适的范围内，通过建立数学模型解决以下三个问题。

(1)建立只有中心室的一室模型，求解在快速静脉注射、恒速静脉注射（持续时间为τ）和口服或肌肉注射3种给药方式下的血药浓度，并画出血药浓度曲线的图形。

(2)在多次重复快速静脉注射给药方式下，写出血药浓度表达式并作图。

讨论怎样确定时间间隔和给与固定剂量，使血药浓度的变化满足上述要求（实际上为了简化起见，常采用加大首次剂量给药的方式，给出这种情况下的给药方案）。

(3)在多次重复恒速静脉滴注和口服（或肌肉注射）的给药方式下，给出血药浓度的简图，并选择一种方式讨论确定时间间隔和每次给予固定剂量的问题。

二、问题分析
针对问题一，先建立中心室药物含量随时间变化的微分方程模型（[1]），其中药物含量的变化取决于注入速率与中心室自身排泄速率，然后解出药物含量随时间变化的关系式，再根据浓度、体积和药物含量的联系求出血药浓度随时间变化的关系式。

查阅资料得到某种药物的具体实例，运用Matlab软件进行曲线拟合([2])可确定药物动力学参数（[3]），画出血药浓度曲线的图形。

而快速静脉注射、恒速静脉注射和口服或肌肉注射这三种情况是在微分方程不同初值条件下的特解。

对于快速静脉注射，起始血药浓度为药物总含量与中心室体积之比，中心室的给药速率为零；对于恒速静脉注射，起始血药浓度为零，中心室的给药速率为常数，即药物总含量与注射时间τ的比值；对于口服或肌肉注射，起始血药浓度为零，然后以肌肉为另一房室，建立该室药物含量随时间变化的微分方程模型，其中药物含量变化由该房室自身排泄药量组成，代入初值条件即开始注入药物的总含量，求出该房室药物含量随时间变化的关系式。

通过该关系式可求出中心室的给药速率。

把三种不同的初始条件带入微分方程模型中，即可求出三种不同注药方式下血药浓度随时间变化的关系式。

针对问题二，先设定注射药物时间间隔和药物剂量。

采用问题一中一次给药方式下的血药浓度表达式，利用迭代算法算出多次给药方式下血药浓度随时间变化的表达式，用Matlab软件画出血药浓度随时间变化的图像。

通过控制变量法，查阅资料分别赋予给药时间间隔与给药剂量不同的值，运用Matlab软件画出多组血药浓度随时间变化的曲线图，进行分析对比确定给药方案。

针对问题三，在恒速静脉滴注和口服或肌肉注射的多次给药方式下，血药浓度表达式计算方法同问题二。

选择恒速静脉滴注多次给药方式讨论时间间隔和固
定剂量，采用上述相同的控制变量法进行讨论。

三、模型假设
(1)忽略注射到肌肉中的药物转化到中心室所需时间。

(2)假设肌肉注射的药物全部转化到中心室。

(3)假设以肌肉为吸收室的药物排泄速率和转移到中心室的速率为常数。

(4)假设中心室药物排泄速率为常数。

(5)假设中心室的体积在整个过程中保持不变。

(6)假设注入的药物全进入中心室。

四、符号表示
τ恒速静脉注射持续时间
c血药浓度
t时间
v中心室体积
k中心室药物排泄速率
r药物进入中心室速率
k恒速静脉滴注的药物注射速率
1
k肌肉内药物排泄速率
2
k肌肉到中心室的转移速率
3
t恒速滴注时的时间间隔
m初始注入药物剂量
()t x中心室内t时刻的药物含量
()t y肌肉内t时刻的药物含量
五、模型的建立与求解
本文在只有中心室的条件下，考虑中心室内血药浓度随时间的变化，建立微分方程模型。

5.1模型的建立
按照快速静脉注射、恒速静脉滴注和口服或肌肉注射3种给药方式，不同的初值对应微分方程不同的解，分别得到3种给药方式下血药浓度随时间变化的表达式。

5.1.1针对问题一
考虑中心室内药物含量随时间变化的关系，建立两者之间关系的微分方程。

快速静脉注射、恒速静脉滴注、口服或肌肉注射这三种不同的给药方式是微分方
程在不同初值条件下的特解，其中药物含量的变化由该时刻注入的药量和中心室排除的药量组成，变化过程如图1所示：
图 1 一室模型示意图
中心室内药物含量随时间变化的微分方程为： ()()t kx r dt
t dx -= (1) 运用常微分方程相关知识解出中心室内药物含量与时间变化的关系式，又由于血药浓度等于药物含量与中心室体积之比，即：
()()v
t x t c = (2) 由(1)和(2)可计算出中心室内血药浓度随时间变化的关系式。

当快速静脉注射时，起始血药浓度v
m c 0=，起始的注射速率为0； 当恒速静脉滴注时，起始血药浓度0=c ，起始的注射速率为1k ；
当口服或肌肉注射时，起始血药浓度0=c 。

假定以肌肉为房室，建立该房室内药物含量随时间变化的微分方程，其中某时刻的药物含量变化由该时刻排除药物量所决定，变化过程如图2所示：
图 2 药物经吸收室进入中心室示意图
肌肉内药物含量随时间变化的微分方程为：
()()()⎪⎩⎪⎨⎧=-=0
20m y t y k dt t dy (3)
通过该微分方程可求出肌肉内药物含量随时间变化的关系式，进而算出该种给药方式下肌肉到中心室的转移速率。

最后把不同给药方式下的初始条件带入微分方程的通解中，算出各自给药方式下血药浓度随时间变化的方程式，用Matlab 软件画出不同的血药浓度曲线图
形。

快速静脉注射的多次重复给药方式下,在上述一次给药模型的基础上利用迭代得出表达式，从而确定给药时间间隔和剂量。

5.1.2针对问题二
设定快速静脉注射的时间间隔为0t ，注射的剂量为0m 。

采用问题一中一次给药方式下的血药浓度表达式进行迭代，求出多次重复给药方式下的血药浓度表达式，运用Matlab 软件画出其图形。

通过控制变量法，确定适当地时间间隔和固定剂量。

首先固定给药剂量，根据经验和查阅资料取6组不同的时间间隔，分别用Matlab 软件画出各自血药浓度随时间变化的曲线图，对比分析确定出最佳时间间隔。

其次选取适当时间间隔，取3组不同的给药剂量，分别用Matlab 软件画出各自血药浓度随时间变化的曲线图，对比分析可找出最佳给药剂量。

以下选择恒速滴注进行讨论，相对以上快速注射，其应根据给药持续时间的不同分段求解。

5.1.3针对问题三
在恒速滴注和口服或肌肉注射多次重复给药方式条件下，求解两者通式的方法与问题二相同，并运用Matlab 软件画出各自给药方式下血药浓度曲线图。

选择恒速静脉滴注的给药方式确定其时间间隔和固定剂量，其确定方法与上述控制变量法相同。

5.2模型的求解
本文建立的是含参数常微分模型，首先应进行参数的确定。

5.2.1药物动力学参数的确定
以某种药为具体事例，查阅相关资料得到相关数据([4]),一室模型快速静脉注射下的血药浓度数据如下表1所示：(t=0注射300mg)
表1 中心室血药浓度—时间数据表
t(h) 0.25 0.5 1 1.5 2 3 4
6 8 c(ug/ml) 19.21 18.15 15.36 14.10 12.89 9.32
7.45 5.24 3.01
用Matlab 软件进行最小二乘曲线拟合（程序代码见附录1），得
中心室体积v =15022ml ，中心室排泄率k =0.2347，拟合图像，如图3所示：
图 3 药物浓度随时间变化的曲线图
在本文中设常数2,1,221===k k τ
对于上述所列微分方程，利用常数变易法进行求解。

5.2.2针对问题一
运用微分方程相关知识，求得微分方程(1)的通解为：
()kt e c k r t x -+=1 (4) 由公式(2)和公式(3)可求得中心室内血药浓度随时间变化的关系式为：
()⎪⎭
⎫ ⎝⎛+=-kt e c k r v t c 1 (5) 当快速静脉注射时，起始血药浓度v
m c 0=
，注射速率01=k ，初值条件代入公式(5)中求得该种给药方式下血药浓度随时间变化的表达式为： ()v
e m t c kt
-=0 (6) 运用软件画出其曲线图形（程序编码见附录2），如图4所示：
图 4 快速注射下血药浓度随时间变化图
当恒速静脉滴注时，起始血药浓度0=c ，注射速率τ0
1m k =，初值条件代入
公式(5)中求出该种给药方式下血药浓度随时间变化的表达式，当药物作用时间τ>t 时，表达式为：
()()()⎪⎪⎩
⎪⎪⎨⎧-+=-=----)(1][10kt kt kt t k e k r e m k r v t c e k r e k r v t c τ t
t ≤≤≤ττ0 (7)
运用Matlab 软件画出其曲线图形（程序代码见附录3），如图5所示：
图 5 恒速滴注下血药浓度随时间变化图
当口服或肌肉滴注时，起始浓度0=c ，由微分方程(5)求出肌肉内药物含量随时间变化的表达式为：
()t k e m t y 20-=
肌肉到中心室的转移速率为：
t k e m k k 2023-=
初值条件代入(5)中求出该种给药方式下血药浓度随时间变化的表达式为：
()()
kt t k e e k k k m t c ----=2220 运用Matlab 软件画出其曲线图形（程序代码见附录4），如图6所示：
图 6 口服或肌肉注射下血药浓度随时间变化图
由问题一解出的表达式，利用迭代算法得出快速注射多次给药方式下的表达式。

5.2.3针对问题二
第一个时间间隔0t 内血液浓度随时间变化的表达式（由公式(6)得出）为：
()v e m t c kt
-=0 利用迭代方法，得出
第二个时间间隔0t 内血液浓度随时间变化的表达式为：
()()
()0010t t k kt e e m v t c ---+= 第三个时间间隔0t 内血液浓度随时间变化的表达式为：
()()000220]1)[(t t k kt kt e e e m v t c ----++=
第四个时间间隔0t 内血液浓度随时间变化的表达式为：
()()00003230]1)()[(t t k kt kt kt e e e e m v t c -----+++=
根据前四个周期内血液浓度随时间变化的表达式归纳总结得出血液浓度随时间变化的表达式为：
()()[]}11{10
010kt nkt t n t k e e e m v t c ------= () 4,3,2=n 运用Matlab 软件画出其血液浓度随时间变化的曲线图（程序代码见附录5），如图7所示：
图 7 快速注射多次重复给药方式下血药浓度随时间变化图
在讨论时间间隔和剂量时，设定最佳血药浓度变化范围，固定剂量0m 的大小，用Matlab 软件画出在最佳血药浓度范围内不同时间间隔0t 下的图像。

当20=t 时，血药浓度随时间变化曲线图如图8所示：
当40=t 时，血药浓度随时间变化曲线图如图9所示：
当60=t 时，血药浓度随时间变化曲线图如图10所示：
当80=t 时，血药浓度随时间变化曲线图如图11所示：
当100=t 时，血药浓度随时间变化曲线图如图12所示：
当160=t 时，血药浓度随时间变化曲线图如图13所示：
图8周期为2时血药浓度变化 图9周期为4时血药浓度变化
图10周期为6时血药浓度变化 图11周期为8时血药浓度变化
图12 周期为10时血药浓度变化 图13 周期为16时血药浓度变化
由图8,9,10,11,12,13分析知，血药浓度分布在最佳浓度范围内的曲线越密集，此时的时间间隔越佳。

固定时间间隔0t ，改变剂量0m 的大小，用Matlab 软件画出血药浓度在不同剂量0m 下的曲线图像。

当mg m 2000=时，血药浓度随时间变化曲线图如图14所示：
当mg m 3000=时，血药浓度随时间变化曲线图如图7所示：
当mg m 4000=时，血药浓度随时间变化曲线图如图15所示：
图14剂量为200时血药浓度变化 图15剂量为400时血药浓度变化
由图14,7,15分析知，血药浓度分布在最佳浓度范围内的曲线越密集，此时注射剂量越佳。

由问题一解出的表达式，利用迭代算法得出恒速滴注和口服或肌肉注射多次给药方式下的表达式。

5.2.4针对问题三
对于恒速滴注多次重复给药方式，采用迭代方法。

第一个时间间隔内血药浓度随时间变化的表达式为：
()()()⎪⎪⎩⎪⎪⎨⎧-=-=---][1)(1kt t k kt e k
r e k r v t c e k r k r v t c τ 00t t t ≤≤≤≤ττ 第二个时间间隔内血药浓度随时间变化的表达式为：
()()()()()()()⎪⎪⎩⎪⎪⎨⎧-+-=-+-=--+-------}]{[1}]{[10000000t t k t t k kt t k t t k kt t k e k
r e k r e k r e k r v t c e k r k r e k r e k r v t c τττ 00222t t t t ≤≤≤≤ττ 第三个时间间隔内血药浓度随时间变化的表达式为：
()()()()()()()⎪⎪⎩⎪⎪⎨⎧-+-=-+-=--+-------}]{[1}]{[10000000222t t k t t k kt t k t t k kt t k e k
r e k r e k r e k r v t c e k r k r e k r e k r v t c τττ003332t t t t ≤≤≤≤ττ 第n 个时间间隔内血药浓度随时间变化的表达式为：
()()()()()()()()()()⎪⎪⎩
⎪⎪⎨⎧-+--=-+--=----+--------})])(1{[(1})])(1{[(10000000111t n t k t n t k kt t k t n t k kt t k e k r e k r e k r e k r n v t c e k r k r e k r e k r n v t c τττ 00t t t ≤≤≤≤ττ
运用Matlab 软件画出恒速滴注多次重复给药方式下血药浓度随时间变化曲线图（程序代码见附录6），如图16所示：
图 16 恒速滴注多次重复给药方式下血药浓度随时间变化图
对于口服或肌肉注射多次重复给药方式，采用迭代方法。

可得出：
第一个时间间隔内血药浓度随时间变化的表达式为：
())]([122
20kt t k e e k k k m v t c ----= 00t t ≤≤ 第二个时间间隔内血药浓度随时间变化的表达式为：
()())]()([100022
20220t t k kt t k e k k k m e e k k k m v t c -----+--= 002t t t ≤≤ 第三个时间间隔内血药浓度随时间变化的表达式为：
()())]()(2[100022
20220t t k kt t k e k k k m e e k k k m v t c -----+--= 0032t t t ≤≤ 第n 个时间间隔内血药浓度随时间变化的表达式为：
()())]()()1[(100022
20220t t k kt t k e k k k m e e k k k m n v t c -----+---= () 4,3,2=n 运用Matlab 软件画出口服或肌肉注射多次重复给药方式下血药浓度随时间
变化曲线图（程序代码见附录7），如图17所示：
图 17 口服或肌肉注射多次重复给药方式下血药浓度随时间变化图
本文讨论在恒速滴注给药方式下的时间间隔和固定剂量。

在讨论时间间隔和剂量时，设定最佳血药浓度变化范围，固定剂量0m 的大小，用Matlab 软件画出在最佳血药浓度范围内不同时间间隔0t 下的图像。

当40=t 时，血药浓度随时间变化的曲线图如图18所示：
当60=t 时，血药浓度随时间变化的曲线图如图16所示：
当80=t 时，血药浓度随时间变化的曲线图如图19所示：
图 18 周期为4时血液浓度变化 图19 周期为8时血液浓度变化 对图18，16，19进行分析可知，血药浓度分布在最佳浓度范围内的曲线越密集，此时的时间间隔越佳。

由表达式知血药浓度与剂量无关。

对以上求解过程及结果对比分析，可得不同给药方式下血药浓度变化规律。

5.3结论分析
本模型讨论了不同给药方式下血药浓度随时间变化的规律，给出了不同给药方式下，确定药物给药时间间隔和剂量的方法，对确立合理的给药方案有很大的参考价值，为临床医学提供了理论依据。

六、模型的评价与推广
本文忽略了药物在进入中心室过程中，体积、排泄速率等因素的变化 ,建立的是简化的线性微分方程模型。

6.1模型的评价
优点：
(1)该模型简单易懂，便于Matlab编程计算。

(2)在确定时间间隔和剂量时，进行多组数据对比分析。

缺点：
(1)建立模型的过程中，简化了中心室对药物吸收、排除的过程，导致结论与实际有一定的差距。

(2)在模型的求解过程中，只对时间间隔和固定剂量进行定性分析，没有给出确定的数值。

6.2模型的推广
(1)在实际生活中为了简化起见，常采用加大首次剂量给药的方式。

假设第一次给药量达到了机体所能承受的最大剂量（即此时对应机体所能承受的最大血药浓度），第二次开始注入恒定的剂量，且其快速注入后恰能达到机体所能承受的最大剂量。

此时血药浓度随时间变化的函数是一个以时间间隔为周期的周期函数，其中每个周期内血药浓度随时间变化的表达式与问题一在该种给药方式下血药浓度随时间变化表达式一样，可运用Matlab软件画出此种给药方式下血药浓度曲线图。

由于时间关系，在此不做详细讨论。

(2)在实际情况中，许多因素是离散变化的而不是均匀连续的，如机体在一天中的不同时刻，体内的新陈代谢速率不同，对药物吸收和排除的速率也不同。

可将本文中的线性模型推广到非线性模型,如常见的非线性模型
()
c
k
c
k
dt
t
dc
+
-
=
4
5。

(3)在生活应用上，可以将该模型应用在植物中残留农药的测定、血液中酒精浓度的测定，其对临床医学的研究有很大的参考价值。

参考文献
[1] 唐焕文，数学建模引论，北京：高等教育出版社，2001.
[2] 张智丰，数学软件与大学数学实验，北京：高等教育出版社，2013.
[3] 王广基，药物代谢动力学，北京：化学工业出版社，2005.
[4] 刘建平，生物药剂学与药物动力学，北京：人民卫生出版社，2011.
附录
附录1
曲线拟合程序：
>> t=[0.25 0.5 1 1.5 2 3 4 6 8];
>> c=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01];
>> y=log(c);
>> a=polyfit(t,y,1)
a =
-0.2347 2.9943
>> k=-a(1)
k =
0.2347
>> v=exp(log(300000)-a(2))
v =
1.5022e+04
z=polyval(a,t);
plot(t,c,'k+',t,z,'r')
附录2
一次快速注射程序：
>> t=[0:1/100:50];
>> c=300000./((1.5022e+04).*exp(0.2347.*t));
plot(t,c)
附录3
一次恒注程序：
function c=funx(t)
if 0<=t&&t<2
c=(1/(0.2347*15022)).*(1-exp(-0.2347.*t));
elseif t>=2
c=(1/(0.2347.*15022)).*(exp(0.2347.*(2-t))-exp(-0.2347.*t));
end
调用函数：
t=0:0.1:30;
for i=1:length(t)
c(i)=fun1(t(i));
end
plot(t,c)
附录4
一次口服或肌肉注射的图像：
t=0:0.1:24;
>> c=(300000*1./(15022*(1-0.2347))).*(exp(-0.2347.*t)-exp(-1.*t)); >> plot(t,c)
附录5
快速注射多次周期为6的程序：
t=0:0.1:6;
c1=(300000/15022).*exp(-0.2347.*t);
plot(t,c1)
hold on
t=6:0.1:12;
c2=(300000/15022)*(exp(-0.2347*4)+1).*exp(-0.2347.*(t-6));
plot(t,c2)
hold on
t=12:0.1:18;
c3=(300000/15022)*((exp(-0.2347*4))^2+exp(-0.2347*4)+1).*exp(-0.2347. *(t-12));
plot(t,c3)
hold on
t=18:0.1:24;
c4=(300000/15022)*((exp(-0.2347*4))^3+(exp(-0.2347*4))^2+exp(-0.2347* 4)+1).*exp(-0.2347*(t-18));
plot(t,c4)
hold on
t=24:0.1:30;
c5=(300000/15022)*((exp(-0.2347*4))^4+(exp(-0.2347*4))^3+exp(-0.2347* 4)^2+exp(-0.2347*4)+1).*exp(-0.2347*(t-24));
plot(t,c5)
附录6
恒速注射多次周期为6的程序
t=0:0.1:2;
c1=(1/15022*0.2347).*(1-exp(-0.2347.*t)) ;
plot(t,c1)
hold on
t=2:0.1:6;
c2=(1/15022*0.2347)*(1-exp(-0.2347*2)).*exp(-0.2347.*(t-2));
plot(t,c2)
hold on
t=6:0.1:8;
c3=(1/15022*0.2347).*((1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp(-0 .2347.*(t-6))));
plot(t,c3)
hold on
t=8:0.1:12;
c4=(1/15022*0.2347).*((1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp(-0 .2347*2)).*exp(-0.2347.*(t-6-2)));
plot(t,c4)
hold on
t=12:0.1:14;
c5=(1/15022*0.2347).*(2*(1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp( -0.2347.*(t-12))));
plot(t,c5)
hold on
t=14:0.1:18;
c6=(1/15022*0.2347).*(2*(1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp( -0.2347*2)).*exp(-0.2347.*(t-12-2)));
plot(t,c6)
hold on
t=18:0.1:20;
c7=(1/15022*0.2347).*(3*(1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp( -0.2347.*(t-18))));
plot(t,c7)
hold on
t=20:0.1:24;
c8=(1/15022*0.2347).*(3*(1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp( -0.2347*2)).*exp(-0.2347.*(t-18-2)));
plot(t,c8)
hold on
t=24:0.1:26;
c9=(1/15022*0.2347).*(4*(1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp( -0.2347.*(t-24))));
plot(t,c9)
hold on
t=26:0.1:30;
c10=(1/15022*0.2347).*(4*(1-exp(-0.2347*2))*exp(-0.2347*(6-2))+(1-exp (-0.2347*2)).*exp(-0.2347.*(t-24-2)));
plot(t,c10)
附录7
口服或肌肉注射多次周期为4的编程：
t=0:0.1:4;
c1=(300000*1./(15022*(1-0.2347))).*(exp(-0.2347.*t)-exp(-1.*t)); plot(t,c1)
hold on
t=4:0.1:8;
c2=(300000*1./(15022*(1-0.2347))).*((exp(-0.2347*4)-exp(-1*4))+(exp(-0.2347.*(t-4))-exp(-1.*(t-4))));
plot(t,c2)
hold on
t=8:0.1:12;
c3=(300000*1./(15022*(1-0.2347))).*(2*(exp(-0.2347*4)-exp(-1*4))+(exp (-0.2347.*(t-8))-exp(-1.*(t-8))));
plot(t,c3)
hold on
t=12:0.1:16;
c4=(300000*1./(15022*(1-0.2347))).*(3*(exp(-0.2347*4)-exp(-1*4))+(exp (-0.2347.*(t-12))-exp(-1.*(t-12))));
plot(t,c4)
hold on
t=16:0.1:20;
c5=(300000*1./(15022*(1-0.2347))).*(4*(exp(-0.2347*4)-exp(-1*4))+(exp (-0.2347.*(t-16))-exp(-1.*(t-16))));
plot(t,c5)
hold on。