function exmr1r2b
  c10 = 10; c20 = 31.5; V0 = 0.100;
  
  param.k = 4.79e-6*3600;
  param.A = 1;
  param.clim = 319;
  param.c1f = 100;
  param.c2f = 10;
  param.R1 = 0.80;
  param.R2 = 0.20;
  
  alpha0 = 0;
  alpha1 = 1; %%% CVD control
  alpha2 = 0;
  
%  St = @(cAclim,r1) log((1/(cAclim*r1)-(1-r1)/r1))-1/(1-(1-r1)*cAclim);
%  param.copt = fzero(@(cAclim) St(cAclim,param.R1),1/exp(1)) * param.clim;
  param.copt = param.clim / exp(1); %%% CVD step
  
  AbsTol = 1e-13;
  RelTol = 1e-12;
  options1 = odeset('Events',@event1,'AbsTol',AbsTol,'RelTol',RelTol);
  state01      = [c10;c20;V0;0];
  tspan1   = [0 10000];
  [time1,state1,tt1] = ode15s(@model,tspan1,state01,options1,alpha0,param);
  
  options2 = odeset('Events',@event2,'AbsTol',AbsTol,'RelTol',RelTol);
  state02      = state1(end,:);
  tspan2   = [time1(end) 10000];
  [time2,state2,tt2] = ode15s(@model,tspan2,state02,options2,alpha1,param);
  
  options3 = odeset('Events',@event3,'AbsTol',AbsTol,'RelTol',RelTol);
  state03      = state2(end,:);
  tspan3   = [time2(end) 10000];
  [time3,state3,tt3] = ode15s(@model,tspan3,state03,options3,alpha2,param);

  state = [state1;state2;state3];
  time = [time1;time2;time3];
  alpha = [alpha0*ones(length(time1),1); alpha1*ones(length(time2),1); alpha2*ones(length(time3),1)];
  figure
  subplot(2,2,1)
  plot(state(:,1),state(:,2))
  hold on
  plot(c10,c20,'go',param.c1f,param.c2f,'rx')
  xlabel('c1-unit')
  ylabel('c2-unit')
  
  subplot(2,2,2)
  plot(time, 1000*state(:,3))
  xlabel('time')
  ylabel('V-unit')
  
  subplot(2,2,3)
  plot(time, 1000*state(:,4))
  xlabel('time')
  ylabel('Vw-unit')
  x = axis; axis([x(1), x(2), -0.2, x(4)]);
  
  subplot(2,2,4)
  plot(time, alpha)
  xlabel('time')
  ylabel('alpha')
  x = axis; axis([x(1), x(2), -0.05, 1.10]);
  str = ['Final time is ' num2str(time(end)) ' h'];
  disp(str)

  set(gcf, 'PaperUnits', 'centimeters') 
  set(gcf, 'PaperPosition', [3 10 11 8.25])
  print -depsc 'exmr1r2b.eps'

function fun = model(t,x,u,param)
  c1 = x(1); c2 = x(2); V = x(3); alpha = u;
  R1 = param.R1; R2 = param.R2; 
  
  q = param.k*param.A*log(param.clim/(c1*R1) - (1-R1)/R1);
  
  dc1dt = c1*q*(R1 - alpha)/V;
  dc2dt = c2*q*(R2 - alpha)/V;
  dVdt  = (alpha-1)*q;
  daqdt = alpha*q;
  
  fun = [dc1dt;dc2dt;dVdt;daqdt];
  
function [value,isterminal,direction] = event1(t,y,u,param)
  
  value      = y(1) - param.copt;
  isterminal = 1;
  direction  = 1;
  
function [value,isterminal,direction] = event2(t,y,u,param)
  
  value      = param.c1f/y(1) - (param.c2f/y(2))^(param.R1/param.R2);
  isterminal = 1;
  direction  = 0;
  
  function [value,isterminal,direction] = event3(t,y,u,param)
  
  value      = y(1) - param.c1f;
  isterminal = 1;
  direction  = 0;