function [F,Pcal]=fun_main(Do,k,t_save,t_value,flag)
%此程序返回目标函数的值.
%定义计算区域尺寸
LW=5;
LO=10-5;
NN=length(t_save);
%定义网格
Nw=100;
No=100;
dt=100;
Tmax=15*3600;
dzw=LW/Nw;
dzo=LO/No;
Wgrid=LW/2:-dzw:-LW/2;
Ogrid=LO/2:-dzo:-LO/2;
%构造扩散项差分矩阵
Dz2W_unify=Dif2c(Nw);Dz2W_unify=sparse(Dz2W_unify);Dz2W=Dz2W_unify/((LW/2)^2);
Dz2O_unify=Dif2c(No);Dz2O_unify=sparse(Dz2O_unify);Dz2O=Dz2O_unify/((LO/2)^2);
[Dz1W_unify,~]=Dif(Nw);Dz1W_unify=sparse(Dz1W_unify);Dz1W=Dz1W_unify/(LW/2);
[Dz1O_unify,~]=Dif(No);Dz1O_unify=sparse(Dz1O_unify);Dz1O=Dz1O_unify/(LO/2);
Iw=speye(Nw+1);
Io=speye(No+1);
%初始条件
Dw=1e-5;
% Do=4.5e-5;
% k=3;
Ci=1.05e-3;
Cw0=Ci*ones(Nw+1,1);
Co0=zeros(No+1,1);
CWold=Cw0;
COold=Co0;
V0=100;m0=71.93;Mco2=44;
%求解扩散问题Dc/dt= Dw*Dz2W*c Dc/dt= Dz2O*c
%隐式差分格式 (cn+1-cn)/dt= Dw*Dz2W*cn+1
%cn= (1-Dw*Dz2W*dt)*cn+1
A1=(Iw(2:end-1,:)-Dw*Dz2W*dt);
A2=(Io(2:end-1,:)-Do*Dz2O*dt);
Ai=[A1,zeros(Nw-1,No+1);zeros(No-1,Nw+1),A2];
Ab=[Dz1W(end,:), zeros(1,No+1);...
zeros(1,Nw+1), Dz1O(1,:);...
-k*Iw(1,:), Io(end,:);...
Dw*Dz1W(1,:),-Do* Dz1O(end,:)];
AA=[Ai;Ab];
% t_save=[0,5*3600,10*3600,15*3600];
mycolor=jet(length(t_save));
ii=1;
for t=0:dt:Tmax
BB=[CWold(2:end-1);COold(2:end-1);zeros(4,1)];
Cnew=AA\BB;
CWnew=Cnew(1:Nw+1);
COnew=Cnew(Nw+2:Nw+No+2);
CWold=CWnew;
COold=COnew;
if length(find(t_save==t))>0&flag==1
figure(3)
plot(LW/2+Wgrid,CWnew,'Color',mycolor(ii,:))
hold on;
plot(LO/2+LW+Ogrid,COnew,'Color',mycolor(ii,:))
xlabel('x');
ylabel('C');
nco2=-trapz(Ogrid,COnew);
Cave=nco2/V0;
rou=(Mco2*nco2+m0)/V0;
Pcal(ii)=(rou-0.6874-3.1*Cave)/7.97e-3;
ii=ii+1;
elseif length(find(t_save==t))>0&flag==0
nco2=-trapz(Ogrid,COnew);
Cave=nco2/V0;
rou=(Mco2*nco2+m0)/V0;
Pcal(ii)=(rou-0.6874-3.1*Cave)/7.98e-3;
ii=ii+1;
end
end
F=(1/NN)*norm(Pcal-t_value);
end