function example_dirac_secondorder_vectmult()

COMPILATION

    [exdir,~,~]=fileparts(which('example_dirac_secondorder_vectmult.m'));
    % compile the model
    amiwrap('model_dirac_secondorder_vectmult','model_dirac_secondorder_vectmult_syms',exdir,2)
Generating model struct ...
x | k | p | deltax | xdot | deltaxdot | ddeltaxdx | ddeltaxdp | ddeltaxdt | root | drootdx | sx | drootdp | drootdt | dtaudp | sroot | stau | deltasx | sigma_y | dsigma_ydp | y | dydp | sigma_z | dsigma_zdp | z | dzdp | Parsing model struct ...
z | Generating C code ...
deltasx | deltax | dsigma_ydp | dsigma_zdp | dydp | dzdp | root | sigma_y | sigma_z | stau | xdot | y | z | headers | wrapfunctions | headers | wrapfunctions | Compiling mex file ...
amici | Building with 'Xcode with Clang'.
MEX completed successfully.
Building with 'Xcode with Clang'.
MEX completed successfully.
amici | Building with 'Xcode with Clang'.
MEX completed successfully.
Building with 'Xcode with Clang'.
MEX completed successfully.

SIMULATION

    % time vector
    t = linspace(0,3,1001);
    p = [1;0.5;2;3];
    k = [];
    v = [0.7;0.3;1.4;0.1];

    options = amioption('sensi',0,...
        'maxsteps',1e4);

    % load mex into memory
    [msg] = which('model_dirac_secondorder_vectmult'); % fix for inaccessability problems
    options.sensi = 2;
    sol = simulate_model_dirac_secondorder_vectmult(t,log10(p),k,[],options,v);

FORWARD SENSITIVITY ANALYSIS

    options.sensi = 2;

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

FINITE DIFFERENCES

    options.sensi = 1;

    eps = 1e-4;
    xi = log10(p);
    for ip = 1:4;
        xip = xi;
        xip(ip) = xip(ip) + eps;
        solp = simulate_model_dirac_secondorder_vectmult(t,xip,k,[],options);
        s2x_fd(:,:,ip) = sum(bsxfun(@times,(solp.sx - sol.sx)/eps,permute(v,[3,2,1])),3);
        s2y_fd(:,:,ip) = sum(bsxfun(@times,(solp.sy - sol.sy)/eps,permute(v,[3,2,1])),3);
    end

PLOTTING

    figure
    c_x = get(gca,'ColorOrder');
    for ip = 1:4
        subplot(4,2,ip*2-1)
        hold on
        for ix = 1:size(sol.x,2)
            plot(t,sol.s2x(:,ix,ip),'.-','Color',c_x(ix,:))
            plot(t,s2x_fd(:,ix,ip),'--','Color',c_x(ix,:))
        end
        ylim([-10,10])
        legend('x1','x1_{fd}','x2','x2_{fd}','Location','NorthEastOutside')
        legend boxoff
        title(['state sensitivity for p' num2str(ip)])
        xlabel('time t')
        ylabel('x')
        box on

        subplot(4,2,ip*2)
        plot(t,abs(sol.s2x(:,:,ip)-s2x_fd(:,:,ip)),'r--')
        legend('error x1','error x2','Location','NorthEastOutside')
        legend boxoff
        title(['state sensitivity for p' num2str(ip)])
        xlabel('time t')
        ylabel('error')
        ylim([1e-12,1e0])
        set(gca,'YScale','log')
        box on
    end
    set(gcf,'Position',[100 300 1200 500])

    figure
    for ip = 1:4
        subplot(4,2,ip*2-1)
        hold on
        for iy = 1:size(sol.y,2)
            plot(t,sol.s2y(:,iy,ip),'.-','Color',c_x(iy,:))
            plot(t,s2y_fd(:,iy,ip),'--','Color',c_x(iy,:))
        end
        ylim([-10,10])
        legend('y1','y1_{fd}','Location','NorthEastOutside')
        legend boxoff
        title(['observable sensitivity for p' num2str(ip)])
        xlabel('time t')
        ylabel('y')
        box on

        subplot(4,2,ip*2)
        plot(t,abs(sol.s2y(:,:,ip)-s2y_fd(:,:,ip)),'r--')
        legend('error y1','Location','NorthEastOutside')
        legend boxoff
        title(['observable sensitivity for p' num2str(ip)])
        xlabel('time t')
        ylabel('error')
        ylim([1e-12,1e0])
        set(gca,'YScale','log')
        box on
    end
    set(gcf,'Position',[100 300 1200 500])

    drawnow
end