function [dirk] = get_DIRK_coefs(torder,stages)

if nargin < 2
    stages = torder - 1;
end

switch torder
    case 1
        dirk.lambda = 1;
        dirk.aij = 1;
        dirk.dti = 1;
        dirk.bi = 1;
    case 2
        alpha = 1 - 1 / 2 * sqrt(2);
        dirk.lambda = alpha;
        dirk.aij = [alpha, 0;
            1-alpha, alpha];
        dirk.dti = [alpha, 1];
        dirk.bi = dirk.aij(2,:);
        
    case 3
        if stages == 2
            alpha = 1/2/sqrt(3);
            dirk.lambda = 1/2 + alpha;
            dirk.aij = [dirk.lambda, 0;
                -1/sqrt(3), dirk.lambda];
            dirk.dti = [1/2 + alpha, 1/2 - alpha];
            dirk.bi = [1/2,1/2];
        else
            if stages ~=3
                display('Using 3 stages for 3rd order accuracy')
            end
            % alpha = roots([1,-3,3/2,-1/6]);
            % alpha = alpha(logical((alpha < 1/2) .* (alpha > 1/6)))
            alpha = 0.435866521508459;
            
            dirk.lambda = alpha;
            dirk.aij = [alpha, 0, 0;
                (1 + alpha) / 2 - alpha, alpha , 0;
                -(6 * alpha^2 - 16 * alpha + 1) / 4, ... cont
                (6 * alpha^2 - 20 * alpha + 5) / 4, alpha];
            dirk.dti = [alpha, (1 + alpha) / 2, 1];
            dirk.bi = dirk.aij(3,:);
        end
        
    case 4
        alpha = 2 * cos (pi / 18) / sqrt(3);
        dirk.lambda = (1 + alpha) / 2;
        dirk.aij = [dirk.lambda, 0, 0;
            -alpha / 2, dirk.lambda, 0;
            1 + alpha, -(1 + 2 * alpha), dirk.lambda];
        dirk.dti = [(1 + alpha) / 2, 1/2, (1 - alpha / 2)];
        dirk.bi = [1 / (6 * alpha^2), 1 - 1 / (3 * alpha^2),...
            1 / (6 * alpha^2)];
        
        if stages ~= 3
            display('Using 3 stages for 4th order accuracy')
        end
        
    case 5
        a41 = (59765462671 - 2469071899 * sqrt(41)) / 269065110000;
        a42 = (26775007261 + 244199891 * sqrt(41)) / 89286180000;
        a43 = (889326089143 - 203592224167 * sqrt(41)) / 19910818140000;
        dirk.lambda = 1 / 4;
        dirk.aij = [dirk.lambda, 0, 0, 0, 0;
            -1/12, dirk.lambda, 0, 0, 0;
            (73 + 12 * sqrt(41)) / 150, (24 - 19 * sqrt(41)) / 300, ... cont
            dirk.lambda, 0, 0;
            a41, a42, a43, dirk.lambda, 0;
            0, 15/37, (2091 - 879 * sqrt(41)) / 12136, ... cont
            (2091 + 879 * sqrt(41)) / 12136, dirk.lambda];
        dirk.dti = [1/4, 1/6, (49 + sqrt(41))/ 60, (49 - sqrt(41)) / 60, 1];
        dirk.bi = dirk.aij(5,:);
        if stages ~= 5
            display('Using 3 stages for 5th order accuracy')
        end
        
    otherwise
        display('No appropriate DIRK scheme found')
end
            