close all; clear; clc;

%%
figure(1); clf;
y0 = [1 0 0 0 0]';
kAB = 1e3; kBC = 1e3;
kEA = 1; kCE = 10; kDE = kCE;
kCAs = [5e2 5e3 5e4];
for j = 1:3
    kCA = kCAs(j);
    kCDm = @(t) 10*((sin(2*pi*log(t))+1)/2);
    kCd = @(t) kCA+kCE+kCDm(t);  kDd = @(t) kCDm(t)+kDE;
    K = @(t) [-kAB  0   +kCA       0      +kEA ;
        +kAB -kBC  0         0       0   ;
        0   +kBC -kCd(t)   +kCDm(t) 0   ;
        0    0    +kCDm(t) -kDd(t)  0   ;
        0    0    +kCE     +kDE    -kEA];
    
    f = @(t,y) K(t)*y;
    [t,y] = ode15s(f,[1e-5 1e2],y0);
    s = subplot(1,3,j);
    plot(t,y,'-');
    a = gca; a.XScale = 'log';
    grid on;
    xlim([1e-5 1e2]);
    hold on;
    s.XTick = 10.^(-4:2:2);
    if j>1
        s.YTickLabel = [];
    else
        ylabel('concentration')
    end
    t = 10.^linspace(-6,3,1e3);
    l = plot(t,0.925+0.05*((sin(2*pi*log(t))+1)/2),'k-');
    uistack(l,'bottom');
    axis square;
    drawnow
end
