A general BJT model in MATLAB

bjt.m

function [I, J] = bjt(v, Is, N, Bf, Br)
% BJT: a general model for the currents through a BJT and their
% derivatives.
% Ben Holmes, 2017/11/28, GPL 3.0
% Arguments:
%   - v:    [vce; vbe], collector-emitter voltage and base-emitter voltage.
%   - Is:   saturation current
%   - N:    ideality factor
%   - Bf:   forward current gain
%   - Br:   reverse current gain

% Outputs:
%   - I:    [collector current; base current]
%   - J:    jacobian of I with respect to [Vce; Vbe]

% To calculate pnp instead of npn, flip voltage signs and flip current 
% signs.

    % Expand voltages
    vce = v(1);
    vbe = v(2);
    
    % Calculate remaining voltage
    vbc = vbe - vce;
    
    % Pre-multiply N by Vt as they only appear as a product.
    NVt = 25.83e-3*N;
    
    % Simplify reverse current gain factor.
    ar = ((Br + 1)/Br);
    
    % Calculate forward and reverse currents
    If = Is*(exp(vbe/NVt) - 1);
    Ir = Is*(exp(vbc/NVt) - 1);

    % Combine to form collector current and base current
    Ic = If - (Ir*ar);
    Ib = If/Bf + Ir/Br;
    I = [Ic; Ib];

    % Only one derivative is needed, if evaluating vce the polarity is
    % reversed.
    dIf = (If + Is)/NVt;
    dIr = (Ir + Is)/NVt;

    % Calculate the derivatives
    dIcVce = dIr*ar;
    dIcVbe = dIf - (dIr*ar);

    dIbVce = -dIr/Br;
    dIbVbe = dIf/Bf + dIr/Br;

    J = [dIcVce dIcVbe;
        dIbVce dIbVbe];
end