restart; # Define the ODEs for the system odeSys := diff(S(t), t) = z*(1-P-P)+w*C(t)-x*Q*C(t)*S(t)/N-x*Q*Iota(t)*S(t)/N-v*S(t), diff(C(t), t) = z*P+x*Q*C(t)*S(t)/N+x*Q*Iota(t)*S(t)/N-w*C(t)-k*o*C(t)-v*C(t), diff(Iota(t), t) = k*o*C(t)-Iota(t)*v+z*P; collect(add(odeSys[j], j = 1 .. 3), [z, v]); odeTotal := eval(%, Iota(t) = Total(t)-S(t)-C(t)); dsolve({odeTotal, Total(0) = N0}); This formulation comes from copying the 2D version. It is needlessly complicated.Below it I give a more direct way: RHS := evalindets(`~`[rhs]([odeSys]), function, proc (f) options operator, arrow;op(0, f) end proc); RHS:=evalindets(rhs~([odeSys]),function,f->op(0,f)); # The simpler way #the right hand sides E := solve(RHS, {C, Iota, S}); indets(E, RootOf); # Equilibria pol := op([1, 1], %); parnames := indets(RHS, name) minus {C, Iota, S}; parnamesList := convert(parnames, list); ics := S(0) = S0, C(0) = C0, Iota(0) = I0; params := [S0, C0, I0, parnamesList[]]; # Must be given numerical values for an illustration ### Using the parameters option to dsolve/numeric: res := dsolve({ics, odeSys}, numeric, parameters = params); ### Example only. The order is irrelevant if the parameters are given as equations: Example1 := [S0 = 100, C0 = 0, I0 = 0, k = .1, o = 1, v = 0.1e-1, w = .5, x = 1, z = 2, P = .2, P = .3, Q = 1, Q = 2, N = 100]; # Setting parameters res(parameters = Example1); plots:-odeplot(res, [[t, S(t)], [t, C(t)], [t, Iota(t)]], t = 0 .. 500, size = [1200, default], legend = [S, C, Iota]); eval(E, Example1); allvalues(%); # Compare to the result below:` res(500) RHS; # The right hand sides of odeSys s1, c1 := selectremove(has, RHS, S); Again this looks needlessly strange, but is a result of copying 2D math input F := `<,>`(s1, 0); V := `<,>`(eval(c1, z = 0), eval(RHS, z = 0)); The following does the same thing, but looks more natural: F:=<s1,0>; V:=< eval(c1,z=0), eval(RHS,z=0) >; JF1 := VectorCalculus:-Jacobian(F, [C, Iota]); JV := VectorCalculus:-Jacobian(V, [C, Iota]); eval(RHS, {C = 0, Iota = 0, S = N}); JF := eval(JF1, S = N); JF . (-JV)^(-1); ev:=LinearAlgebra:-Eigenvalues(%); Another instance of copying with a complicated version: R0 := `assuming`([max(ev)], [positive]); The better looking way: R0:= max(ev) assuming positive;