System of ODEs Solver

Enter differential equations or matrix coefficients to classify equilibrium points, compute eigenvalues step by step, and generate trajectory phase portraits.

Free to useNo sign-up requiredInstant Results

System of ODEs SolverTry it now — free ▼

Quick examples:

Stable node Center (oscillator) Stable spiral Saddle Repeated (degenerate) Pendulum Van der Pol 3×3 third-order

▦ 2×2 Linear ▦ 3×3 Linear ∿ Nonlinear 2D

Solve \(\dfrac{d}{dt}\!\begin{pmatrix}x\\y\end{pmatrix} = A \begin{pmatrix}x\\y\end{pmatrix}\) where \(A\) is a 2×2 matrix with real entries.

Coefficient matrix A

[

Accepts decimals (e.g. 2.5), negatives (-1), and simple fractions (3/4).

Solve \(\dfrac{d}{dt}\!\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix} = A \begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}\) — get eigenvalues, eigenvectors, and the closed-form general solution.

Coefficient matrix A (3×3)

[

The default is the companion matrix of \( s^3 + 6s^2 + 11s + 6 \), giving eigenvalues \(-1, -2, -3\).

Solve a general 2D nonlinear system numerically: \(x' = f(x, y)\), \(y' = g(x, y)\). Use functions like sin, cos, exp, sqrt, and constants pi, e.

x' = f(x, y)

y' = g(x, y)

No closed-form solution is computed; instead a high-accuracy RK4 integration produces the trajectory and phase portrait.

Initial condition and time span

x(0)

y(0)

time span T

Embed System of ODEs Solver Widget

About System of ODEs Solver

The System of ODEs Solver is an all-in-one differential-equations toolbox for coupled linear and nonlinear systems. Paste a 2×2 or 3×3 coefficient matrix and the tool performs the full eigenvalue / eigenvector analysis, writes the closed-form general and particular solution in LaTeX, classifies the equilibrium at the origin as a saddle, node, spiral, or center, and draws an interactive phase portrait with an animated trajectory. For nonlinear planar systems you can type arbitrary right-hand sides \(f(x,y)\) and \(g(x,y)\) and the tool produces a high-accuracy RK4 phase portrait.

What Is a System of ODEs?

A system of ordinary differential equations couples several unknown functions of a single variable — usually time \(t\) — through their derivatives. In the most compact form,

\[ \mathbf{x}'(t) = \mathbf{F}(t, \mathbf{x}(t)) , \qquad \mathbf{x} \in \mathbb{R}^n \]

When \(\mathbf{F}(t, \mathbf{x}) = A\mathbf{x}\) for a constant matrix \(A\), the system is linear and autonomous — and this is where the theory is at its most beautiful: the entire long-term behavior is determined by the eigenvalues of \(A\).

The Eigenvalue Recipe for Linear Systems

For \(\mathbf{x}' = A\mathbf{x}\) the standard method is:

Compute the characteristic polynomial \(\det(\lambda I - A) = 0\).
Solve for the eigenvalues \(\lambda_1, \lambda_2, \dots\).
For each eigenvalue, find an eigenvector \(v\) by solving \((A - \lambda I) v = 0\).
Assemble the general solution as a linear combination: \(\mathbf{x}(t) = c_1 v_1 e^{\lambda_1 t} + c_2 v_2 e^{\lambda_2 t} + \cdots\).
Fix the constants \(c_i\) by plugging the initial condition \(\mathbf{x}(0)\) into the general solution.

Three Cases for 2×2 Systems

Eigenvalues	General Solution	Portrait
Real distinct \(\lambda_1 \ne \lambda_2\)	\(c_1 v_1 e^{\lambda_1 t} + c_2 v_2 e^{\lambda_2 t}\)	Saddle if signs differ; node otherwise
Complex conjugate \(\alpha \pm i\beta\)	\(e^{\alpha t}[c_1(p\cos\beta t - q\sin\beta t) + c_2(p\sin\beta t + q\cos\beta t)]\)	Spiral (\(\alpha \ne 0\)) or center (\(\alpha = 0\))
Repeated \(\lambda_1 = \lambda_2 = \lambda\)	\(c_1 v e^{\lambda t} + c_2 (tv + w) e^{\lambda t}\)	Degenerate node

The Trace-Determinant Plane

For a 2×2 matrix with trace \(T = a_{11} + a_{22}\) and determinant \(D = a_{11} a_{22} - a_{12} a_{21}\), the whole classification fits in a single diagram:

\[ \begin{array}{l} D < 0 \Rightarrow \text{saddle} \\ D > 0, \; T^2 - 4D > 0 \Rightarrow \text{node (stable if } T < 0\text{)} \\ D > 0, \; T^2 - 4D < 0 \Rightarrow \text{spiral (stable if } T < 0\text{), center if } T = 0 \\ T^2 - 4D = 0 \Rightarrow \text{degenerate node (repeated eigenvalue)} \end{array} \]

This is why the results panel prominently shows \(T\), \(D\), and \(\Delta = T^2 - 4D\) — three numbers are enough to name the equilibrium.

Nonlinear Systems and the Phase Portrait

Most real-world ODEs are nonlinear and have no closed-form solution. The tool handles them by integrating the equations numerically with a 4th-order Runge–Kutta (RK4) method, which has local truncation error \(O(h^5)\) and is the default workhorse for smooth vector fields.

\[ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h/2, y_n + h k_1 / 2) \\ k_3 &= f(t_n + h/2, y_n + h k_2 / 2) \\ k_4 &= f(t_n + h, y_n + h k_3) \\ y_{n+1} &= y_n + \tfrac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} \]

The phase portrait overlays:

A vector field sampled on a 13×13 grid, showing the direction of flow at every point.
The trajectory from your initial condition, drawn in red with an animated orange runner showing the direction of time.
Several seed flow lines from a ring of starting points, giving a global picture of the dynamics.
For 2×2 linear systems, the eigenvector axes (dashed cyan) — these are the invariant directions along which solutions slide exponentially.

How to Use This Solver

Pick a mode — 2×2 Linear, 3×3 Linear, or Nonlinear 2D — via the tabs at the top of the form.
Fill in coefficients or equations. Click any Quick Example to pre-fill a canonical system (stable node, center, saddle, pendulum, Van der Pol, etc.).
Enter an initial condition \((x_0, y_0)\) and a time span \(T\). Typical \(T\) values are 6–20 for oscillators and 3–6 for fast-decaying stable systems.
Click Solve. The full results page appears with the classification, eigenvalues, eigenvectors, closed-form solution (linear modes), animated phase portrait, and time-series plot.
Replay the trajectory using the button under the phase portrait if you'd like to watch the runner traverse the IC curve again.

Worked Example — Damped Harmonic Oscillator

The damped oscillator \(\ddot{x} + 2\zeta \omega \dot{x} + \omega^2 x = 0\) can be rewritten as a 2D system by letting \(y = \dot{x}\):

\[ \begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -\omega^2 & -2\zeta\omega \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \]

For \(\omega = 1\) and \(\zeta = 0.2\) (underdamped), the matrix is \(A = \begin{pmatrix} 0 & 1 \\ -1 & -0.4 \end{pmatrix}\). Trace \(T = -0.4\), determinant \(D = 1\), discriminant \(\Delta = 0.16 - 4 = -3.84 < 0\), so we get a stable spiral with eigenvalues \(-0.2 \pm 0.9798\,i\). The trajectory spirals into the origin, and the time series shows exponentially-decaying sinusoids.

Applications

Mechanical systems — coupled spring-mass systems, pendulums, gyroscopes.
Electrical circuits — RLC networks, op-amp filters, state-space control.
Population dynamics — Lotka–Volterra predator-prey, competing species, epidemiology (SIR, SIS).
Chemical kinetics — reaction networks, Belousov–Zhabotinsky oscillators.
Neuroscience — FitzHugh–Nagumo neuron model, Hodgkin–Huxley reductions.
Control theory — linearized plant models, observer design, stability margins.

Tips & Gotchas

If your trajectory blows up quickly, reduce the time span T — an unstable system can overflow any viewport in a few time units.
For a repeated eigenvalue the solver automatically finds the generalized eigenvector \(w\) by solving \((A - \lambda I)w = v\), so you get the \(tv\) term without manual work.
For nonlinear systems, the vector field arrows also reveal non-origin equilibria as cyan dots — watch the portrait for zero-magnitude regions.
For 3×3 systems there is no phase portrait (3D is hard to display in a 2D page), but the closed-form solution and stability verdict still apply.
Initial conditions and time spans are decoupled from classification: changing them only moves the red trajectory, not the eigenvalue verdict.

Frequently Asked Questions

What is a system of ordinary differential equations?

A system of ordinary differential equations (ODEs) is a set of coupled equations that relate the derivatives of several unknown functions of a single independent variable, typically time. The classic form is \( \mathbf{x}'(t) = F(t, \mathbf{x}(t)) \), where \( \mathbf{x} \) is a vector of states and \(F\) is the vector field. Linear systems can be written compactly as \( \mathbf{x}' = A\mathbf{x} + \mathbf{b} \), and their behavior is determined almost entirely by the eigenvalues of the coefficient matrix \(A\).

How do eigenvalues classify the equilibrium of a 2×2 linear system?

For a 2×2 system \( \mathbf{x}' = A\mathbf{x} \) the origin is classified by the trace \(T\) and determinant \(D\) of \(A\): \(D < 0\) gives a saddle point (unstable); \(D > 0\) with \(T^2 > 4D\) gives a node (stable if \(T < 0\), unstable if \(T > 0\)); \(D > 0\) with \(T^2 < 4D\) gives a spiral (stable if \(T < 0\), unstable if \(T > 0\), a pure center if \(T = 0\)). The borderline \(T^2 = 4D\) produces a degenerate node.

What does the closed-form solution look like when the eigenvalues are complex?

If \(A\) has complex-conjugate eigenvalues \( \alpha \pm i\beta \) with complex eigenvector \( v = p + iq \), the real general solution is \( \mathbf{x}(t) = e^{\alpha t} \left[ c_1 (p \cos\beta t - q \sin\beta t) + c_2 (p \sin\beta t + q \cos\beta t) \right] \). The exponential \(e^{\alpha t}\) controls the amplitude (growing, decaying, or constant), while the sine and cosine handle the rotation.

What happens when the matrix has a repeated eigenvalue?

If the matrix has a repeated eigenvalue \(\lambda\) but only one linearly independent eigenvector \(v\), you also need a generalized eigenvector \(w\) solving \( (A - \lambda I) w = v \). The general solution then takes the form \( \mathbf{x}(t) = c_1 v e^{\lambda t} + c_2 (tv + w) e^{\lambda t} \). If the eigenspace happens to be two-dimensional, the matrix is a scalar multiple of the identity on that invariant subspace and the solution reduces to \( \mathbf{x}(t) = (c_1 v_1 + c_2 v_2) e^{\lambda t} \).

Can this tool solve nonlinear systems symbolically?

The nonlinear mode solves the system numerically using a fourth-order Runge–Kutta (RK4) integrator and plots the phase portrait. Most nonlinear systems have no closed-form solution, so this is the standard approach. You can still read off local behavior near equilibria by linearizing, which the 2×2 linear mode handles — compute the Jacobian at the fixed point and feed it in as \(A\).

What is a phase portrait?

A phase portrait is a geometric picture of the solutions of a 2D system in the \(x\)–\(y\) plane. Each solution traces out a curve called a trajectory, and the collection of trajectories together with the vector field arrows reveals the qualitative behavior: whether solutions spiral in, saddle apart, oscillate, or settle onto equilibria. Phase portraits make the global structure of a system visible at a glance.

System of ODEs Solver

About System of ODEs Solver

What Is a System of ODEs?

The Eigenvalue Recipe for Linear Systems

Three Cases for 2×2 Systems

The Trace-Determinant Plane

Nonlinear Systems and the Phase Portrait

How to Use This Solver

Worked Example — Damped Harmonic Oscillator

Applications

Tips & Gotchas

Frequently Asked Questions

What is a system of ordinary differential equations?

How do eigenvalues classify the equilibrium of a 2×2 linear system?

What does the closed-form solution look like when the eigenvalues are complex?

What happens when the matrix has a repeated eigenvalue?

Can this tool solve nonlinear systems symbolically?

What is a phase portrait?

Further Reading

Calculus:

Math Solver

Numbers & factors

Fractions & division

Powers, roots & logs

Algebra & equations

Calculus

Geometry & trigonometry

Top & Updated:

Download everything as one ZIP with Pro