function example_steadystate

COMPILATION

    [exdir,~,~]=fileparts(which('example_steadystate.m'));
    % compile the model
    amiwrap('model_steadystate','model_steadystate_syms',exdir)
Generating model struct ...
Parsing model struct ...
Generating C code ...
headers | wrapfunctions | Compiling mex file ...
amici | Building with 'Xcode with Clang'.
MEX completed successfully.
Building with 'Xcode with Clang'.
MEX completed successfully.

SIMULATION

    % time vector
    t = linspace(0,300,20);
    p = [1;0.5;0.4;2;0.1];
    k = [0.1,0.4,0.7,1];

    options = amioption('sensi',0,...
        'maxsteps',1e4);
    % load mex into memory
    sol = simulate_model_steadystate(t,log10(p),k,[],options);

    tic
    sol = simulate_model_steadystate(t,log10(p),k,[],options);
    disp(['Time elapsed with cvodes: ' num2str(toc) ])
Time elapsed with cvodes: 0.0040652

ODE15S

    ode_system = @(t,x,p,k) [-2*p(1)*x(1)^2 - p(2)*x(1)*x(2) + 2*p(3)*x(2) + p(4)*x(3) + p(5);
        + p(1)*x(1)^2 - p(2)*x(1)*x(2) - p(3)*x(2) + p(4)*x(3);
        + p(2)*x(1)*x(2) - p(4)*x(3) - k(4)*x(3)];
    options_ode15s = odeset('RelTol',options.rtol,'AbsTol',options.atol,'MaxStep',options.maxsteps);

    tic
    [~, X_ode15s] = ode15s(@(t,x) ode_system(t,x,p,k),t,k(1:3),options_ode15s);
    disp(['Time elapsed with ode15s: ' num2str(toc) ])
Time elapsed with ode15s: 0.10806

PLOTTING

    figure
    c_x = get(gca,'ColorOrder');
    subplot(2,2,1)
    for ix = 1:size(sol.x,2)
        plot(t,sol.x(:,ix),'.-','Color',c_x(ix,:))
        hold on
        plot(t,X_ode15s(:,ix),'d','Color',c_x(ix,:))
    end
    legend('x1','x1_{ode15s}','x2','x2_{ode15s}','x3','x3_{ode15s}','Location','NorthEastOutside')
    legend boxoff
    xlabel('time t')
    ylabel('x')
    box on
    subplot(2,2,2)
    plot(t,abs(sol.x-X_ode15s),'--')
    set(gca,'YScale','log')
    legend('error x1','error x2','error x3','Location','NorthEastOutside')
    legend boxoff
    set(gcf,'Position',[100 300 1200 500])

FORWARD SENSITIVITY ANALYSIS

    options.sensi = 1;
    options.sens_ind = [3,1,2,4];

    sol = simulate_model_steadystate(t,log10(p),k,[],options);

FINITE DIFFERENCES

    eps = 1e-3;

    xi = log10(p);
    for ip = 1:4;
        xip = xi;
        xip(ip) = xip(ip) + eps;
        solp = simulate_model_steadystate(t,xip,k,[],options);
        sx_fd(:,:,ip) = (solp.x - sol.x)/eps;
        sy_fd(:,:,ip) = (solp.y - sol.y)/eps;
    end

PLOTTING

    figure
    for ip = 1:4
        subplot(4,2,ip*2-1)
        hold on
        for ix = 1:size(sol.x,2)
            plot(t,sol.sx(:,ix,ip),'.-','Color',c_x(ix,:))
            plot(t,sx_fd(:,ix,options.sens_ind(ip)),'d','Color',c_x(ix,:))
        end
        legend('x1','x1_{fd}','x2','x2_{fd}','x3','x3_{fd}','Location','NorthEastOutside')
        legend boxoff
        title(['state sensitivity for p' num2str(options.sens_ind(ip))])
        xlabel('time t')
        ylabel('x')
        box on

        subplot(4,2,ip*2)
        plot(t,abs(sol.sx(:,:,ip)-sx_fd(:,:,options.sens_ind(ip))),'--')
        legend('error x1','error x2','error x3','Location','NorthEastOutside')
        legend boxoff
        title(['error of state sensitivity for p' num2str(options.sens_ind(ip))])
        xlabel('time t')
        ylabel('error')
        set(gca,'YScale','log')
        box on
    end
    set(gcf,'Position',[100 300 1200 500])

STEADY STATE SENSITIVITY

    sssens = NaN(size(sol.sx));
    for it = 2:length(t)
        tt = [0,t(it)];
        options.sensi_meth = 'ss';
        solss = simulate_model_steadystate(tt,log10(p),k,[],options);
        sssens(it,:,:) = solss.sx;
        ssxdot(it,:) = solss.xdot;
    end

PLOTTING

    figure
    for ip = 1:4
        subplot(4,2,ip*2-1)
        hold on
        for ix = 1:size(sol.x,2)
            plot(t,sol.sx(:,ix,ip),'.-','Color',c_x(ix,:))
            plot(t,sssens(:,ix,ip),'d-','Color',c_x(ix,:))
        end
        legend('x1','x1_{ss}','x2','x2_{ss}','x3','x3_{ss}','Location','NorthEastOutside')
        legend boxoff
        title(['state steady sensitivity for p' num2str(ip)])
        xlabel('time t')
        ylabel('x')
        box on

        subplot(4,2,ip*2)
        plot(t,abs(sol.sx(:,:,ip)-sssens(:,:,ip)),'--')
        legend('error x1','error x2','error x3','Location','NorthEastOutside')
        legend boxoff
        title(['error of steady state sensitivity for p' num2str(ip)])
        xlabel('time t')
        ylabel('error')
        set(gca,'YScale','log')
        box on
    end
    set(gcf,'Position',[100 300 1200 500])

    figure
    scatter(sqrt(sum((ssxdot./sol.x).^2,2)),sqrt(sum(sum((sol.sx-sssens).^2,2),3)))
    hold on
    plot([1e-15,1e5],[1e-15,1e5],'k:')
    set(gca,'YScale','log')
    set(gca,'XScale','log')
    box on
    axis square
    xlabel('||dxdt/x||_2')
    ylabel('error steady state approximation')
    set(gca,'FontSize',15)
    set(gca,'LineWidth',1.5)
    set(gcf,'Position',[100 300 1200 500])

    drawnow
end