From Equations to Experiments: Modeling ODEs in Julia for Scientific Machine Learning

Over the last few months, I’ve been diving deep into scientific machine learning (SciML), trying to connect the math I’ve learned in class with the kinds of models researchers use to understand real-world systems. I started building and experimenting with ordinary differential equation (ODE) models in Julia, and I quickly realized how powerful — and surprisingly accessible — this workflow can be.

This blog post is my attempt to document that learning process in a way that’s useful to others who are curious about SciML but might be unsure where to start. I’ll walk through how ODEs are represented in Julia, show concrete code examples, and share the intuitions that helped the concepts “click” for me, rather than just restating the theory. Along the way, my goal is to create a resource that feels more like a guided study companion than a polished tutorial, so that other students, self-learners, and early researchers can learn from both my successes and my initial confusion.

If you’re interested in modeling dynamical systems, exploring how differential equations meet machine learning, or just want to see what SciML in Julia looks like in practice, this post is for you.

What is an ODE in practice?

Ordinary Differential Equations are a fundamental class of equations that describe the relationship between functions and their derivatives. It involves one or more functions of one independent variable and it’s derivatives. The term “Ordinary” differentiates these equations from partial differential equations which involve multiple independent variables.

An ODE describes how a state u(t) changes over time according to some rule f:

In Julia, we typically turn this math into a function f!(du, u, p, t) that writes the time derivative into du in-place.

Basic ODE in Julia: exponential decay

Let’s start with a classic example: exponential decay —

Our goal here is to recover the function u(t) that satisfies the differential equation.

How the Julia program helps find u(t)

Mathematically, you can show that the solution is :

In SciML we let the solver approximate this function numerically.

using OrdinaryDiffEq

# du/dt = f(u, p, t)
function decay!(du, u, p, t)
    a = p  # single parameter
    du[1] = -a * u[1]
end

u0    = [1.0]      # initial condition u(0) = 1
p     = 1.5        # decay rate
tspan = (0.0, 5.0) # simulate from t = 0 to 5

prob = ODEProblem(decay!, u0, tspan, p)
sol  = solve(prob, Tsit5())

• ODEProblem bundles the derivative function, initial state, time span, and parameters.
• solve runs a numerical integrator (here Tsit5, a standard adaptive Runge–Kutta method).
• sol(t) lets you query the solution at any time t in the interval.

For solving the problem and deriving the underlying function, we need the following 3 things:

• u0 is the initial condition, i.e., the value of the state u at the starting time. In our case, u0 = [1.0] means we are modeling a quantity that starts at 1 at time t=0.
• p represents the parameters of the model, such as constants that describe the system’s behavior. For exponential decay, p = 1.5 plays the role of the decay rate a, telling the solver how quickly u(t) should decrease over time.
• tspan is the time interval over which we want to solve the ODE, written as a tuple (t_start, t_end). When we set tspan = (0.0, 5.0), we are asking Julia to compute the solution u(t) starting at time 0 and continuing up to time 5.​

How the Julia program helps find u(t)

The OrdinaryDiffEq solver conceptually “walks” along the time axis, repeatedly using decay! to update its estimate of u(t). The result sol is not just a list of numbers—it is an object that behaves like the approximate function t ↦ u(t), so sol(2.0) gives you the numerical value of u(2), and plot(sol) visualizes the whole trajectory.

In other words, instead of manually deriving and coding the closed-form solution, you describe the physical law (the differential equation) and initial condition, and Julia + SciML construct and evaluate the solution function for you.

Inspecting and plotting the solution

using Plots

plot(sol, xlabel = "t", ylabel = "u(t)", title = "Exponential Decay")

plot(sol) automatically uses the time points chosen by the adaptive solver and draws the trajectory.​

When you plot the solution, you see a smooth curve that starts at u(0)=1 and bends downward towards zero without ever actually touching it, which is exactly what a true exponential decay function is.

Conclusion

This post uses the simple exponential decay equation du/dt=−au to show how Julia and the SciML ecosystem turn a mathematical law into an executable model you can simulate, query, and visualize. By defining the differential equation as code, passing it to ODEProblem, and solving it numerically, we recover the familiar decaying curve u(t)=u0*e^−at without ever hard-coding that formula ourselves.