Putting linear ODESystem in matrix-vector form

The symbolic linearization does not yet support analysis points, but this works

eqs = [
    # Ci -> Rih, Rie
    connect(Ci.p, Rie.p),
    connect(Ci.p, Rih.p),

    # heat sources
    # Φh, Φcv -> Ch
    connect(Φhp_c.n, Ch.p),
    connect(Φcv_c.n, Ch.p),

    # Ch -> Rih
    connect(Ch.p, Rih.n),

    # Ce -> Rie, Rea
    connect(Ce.p, Rie.n),
    connect(Ce.p, Rea.p),

    # Ta -> Rea
    connect(Ta.p, Rea.n),
    
    # input to source connections
    connect(Ta_in.output, Ta.V),
    # connect(Φcv_ctrl.output, :Φcv, Φcv_c.I),
    # connect(Φhp_ctrl.output, :Φhp, Φhp_c.I),

    # ground connections
    connect(gnd.g, Ci.n),
    connect(gnd.g, Ce.n),
    connect(gnd.g, Ch.n),
    connect(gnd.g, Φhp_c.p),
    connect(gnd.g, Φcv_c.p),
    connect(gnd.g, Ta.n),
];

components = [
    Ci, Ce, Ch,
    Rih, Rie, Rea,
    Ta_in, Ta,
    Φhp_c, Φcv_c, 
    Φhp_ctrl, Φcv_ctrl, 
    gnd,
];

@named model = ODESystem(eqs, t, systems = components)

ss_matrices, ssys = ModelingToolkit.linearize_symbolic(model, [Φcv_c.I.u, Φhp_c.I.u], [Ci.v, Ce.v, Ch.v])
julia> ss_matrices.A
3×3 Matrix{Num}:
 (-1 / Rih₊R + -1 / Rie₊R) / Ci₊C  1 / (Ci₊C*Rie₊R)                        1 / (Ci₊C*Rih₊R)
                1 / (Ce₊C*Rie₊R)        (-1 / Rea₊R + -1 / Rie₊R) / Ce₊C                  0
      1 / (Ch₊C*Rih₊R)

you may want to further simplify the symbolic entries slightly:

julia> simplify.(ss_matrices.A)
3×3 Matrix{Num}:
             (-Rie₊R - Rih₊R) / (Ci₊C*Rie₊R*Rih₊R)  …   1 / (Ci₊C*Rih₊R)
           1 / (Ce₊C*Rie₊R)                                            0
 1 / (Ch₊C*Rih₊R)