Simultaneous Differential Equations

Basic Concepts

We now describe how to solve two simultaneous first order differential equations of the following form:

y′ = f(x, y, z)          y(x₀) = y₀

z′ = g(x, y, z)          z(x₀) = z₀

Fortunately, the techniques used in Solving Differential Equations using Euler’s Method and Other Methods for Solving Differential Equations are still be applicable. In particular, we can use the following approaches.

Euler’s method

k₁ = hf(x_i, y_i, z_i)           m₁ = hg(x_i, y_i, z_i)

y_i+1 = y_i + k₁          z_i+1 = z_i + m₁

Runge-Kutta (RK2)

k₁ = hf(x_i, y_i, z_i)           m₁ = hg(x_i, y_i, z_i)

k₂ = hf(x_i+h/(2b), y_i+k₁/(2b), z_i+m₁/(2b))           m₂ = hg(x_i+h/(2b), y_i+k₁/(2b), z_i+m₁/(2b))

y_i+1 = y_i + [(1-b)k₁ + bk₂]       z_i+1 = z_i + [(1-b)m₁ + bm₂]

Runge-Kutta (RK4)

k₁ = hf(x_i, y_i, z_i)           m₁ = hg(x_i, y_i, z_i)

k₂ = hf(x_i+h/2, y_i+k₁/2, z_i+m₁/2)           m₂ = hg(x_i+h/2, y_i+k₁/2, z_i+m₁/2)

k₃ = hf(x_i+h/2, y_i+k₂/2, z_i+m₂/2)           m₃ = hg(x_i+h/2, y_i+k₂/2, z_i+m₂/2)

k₄ = hf(x_i+h, y_i+k₃, z_i+m₃)           m₄ = hg(x_i+h, y_i+k₃, z_i+m₃)

y_i+1 = y_i + (k₁+2k₂+2k₃+k₄)/6              z_i+1 = z_i + (m₁+2m₂+2m₃+m₄)/6

Worksheet Function

The Real Statistics Resource Pack provides the following worksheet function that uses Real Statistics lambda capabilities.

DiffEqs(xx, R1, R2, x₀, y₀, z₀, h, ttype, b, Rx, Ry, Rz) = the solution to the differential equations y′ = f(x,y,z) and z′ = g(x,y,z) at the specified x value xx based on the initial values of y(x₀) = y₀ and z(x₀) = z₀

Here R1 contains the expression for f(x,y,z) and R2 contains the expression for g(x,y,z) using the lambda approach described in Differential Equations Support. h is a small positive numeric value as described above. ttype takes the values “ef” for the Euler forward method, “rk” or “rk2” for the Runge-Kutta order 2 method (default ) and “rk4” for the Runge-Kutta order 4 method. If ttype = “rk” or “rk2” then b is as defined above.

The output is a 1 × 2 array with the values y(xx) and z(xx).

If xx – x₀ is not evenly divisible by h, this function uses linear interpolation.

Example 1: Solve the simultaneous differential equations

y′ = 4 cos x – 2 sin x + y – 2z          y(0) = 1

z′ = 5 cos x – 5 sin x – 3y – 4z         z(0) = 2

It turns out that the analytic solution is

y = cos x + sin x          z = 2 cos x

We present the numerical solution in Figure 1 using Euler’s method

Figure 1 – Simultaneous differential equations

Here, range G4:H4 contains the array formula DiffEqs(F4,B$4,B$5,B$6,B$7,B$8,B$9,”ef”). Cell I4 contains the formula =EVALS(B$10,F4) and cell J4 contains =EVALS(B$11,F4).

We can produce a chart of the results as shown in Figure 2.

Figure 2 – Chart of the solution

Using the “rk”, or “rk4” options will provide an even better fit to the actual values.

Examples Workbook

Click here to download the Excel workbook with the examples described on this webpage.

References

Strang, G., Herman, E., Seeburger, P. (2025) Introduction to differential equations
https://math.libretexts.org/Courses/Monroe_Community_College/MTH_211_Calculus_II/Chapter_8%3A_Introduction_to_Differential_Equations/8.1%3A_Basics_of_Differential_Equations

Atkinson, K., Han, W., Stewart, D. (2009) Numerical solutions to ordinary differential equations. Wiley-Interscience
https://homepage.math.uiowa.edu/~atkinson/papers/NAODE_Book.pdf

Hossain, J., Alam, S., Hossain, B. (2017) A study on numerical solutions of second order initial value problems (IVP) for ordinary differential equations with fourth order and Butcher’s fifth order Runge-Kutta methods
http://article.sapub.org/10.5923.j.ajcam.20170705.02.html

Wikipedia (2025) Numerical methods for ordinary differential equations
https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations