衰减方程如下：

http://s16/mw690/003cvAxDzy74YvX8IRV4f&690

matlab程序：
C0=1;
k=0.36;
t=linspace(0,15,90);
dt=t(2)-t(1);
C=C0*exp(-k*t);%解析解
%下面求解数值解
C1=zeros(1,length(t));
C1(1)=C0;
%下面用显式求解
for i=2:length(t)
    C1(i)=C1(i-1)-k*dt*C1(i-1);
end
%下面用隐式求解
C2=C1;
for i=2:length(t)
    C2(i)=C2(i-1)/(1+k*dt);
end
%四阶龙格库塔
lambda=-k;
C3=C1;
for i=2:length(t)
    C3(i)=C3(i-1)*(1+dt*lambda+0.5*(dt*lambda)^2+1/6*(dt*lambda)^3+1/24*(dt*lambda)^4);
end
%半隐半显
C4=C1;
for i=2:length(t)
    C4(i)=C4(i-1)*(1-0.5*k*dt)/(1+0.5*k*dt);
end
plot(t,C,'m--',t,C1,'r--',t,C2,'k--',t,C3,'x',t,C4,'o')
legend('数值解','显式解','隐式解','龙格库塔','半隐半显')
print('-dpng','-r500','C')

fortran程序：

program h1
implicit none
integer::i
real:: C0=1.0,k=0.36,dt=0.1685,lambda=0
real::C(90)=0,C1(90)=0,C2(90)=0,C3(90)=0,C4(90)=0
!解析解
do i=1,90
C(i)=C0*exp(-k*dt*(i-1))
end do
!显式解
C1(1)=C0
do i=2,90
C1(i)=C1(i-1)-k*dt*C1(i-1)
end do
!隐式解
C2(1)=C0
do i=2,90
C2(i)=C2(i-1)/(1+k*dt)
end do
!四阶龙格库塔
C3(1)=C0
lambda=-k
do i=2,90
C3(i)=C3(i-1)*(1+dt*lambda+0.5*(dt*lambda)**2+1.0/6*(dt*lambda)**3+1.0/24*(dt*lambda)**4)
end do
!半隐半显
C4(1)=C0
do i=2,90
C4(i)=C4(i-1)*(1-0.5*k*dt)/(1+0.5*k*dt)
end do
!输出数据
open(1,file='c.txt')
open(2,file='c1.txt')
open(3,file='c2.txt')
open(4,file='c3.txt')
open(5,file='c4.txt')
do i=1,90
write(1,*) C(i)
write(2,*) C1(i)
write(3,*) C2(i)
write(4,*) C3(i)
write(5,*) C4(i)
end do
close(1)
close(2)
close(3)
close(4)
close(5)
end

 http://s2/mw690/003cvAxDzy74YwQ6KM971&690