# ODE returning f(t) = f(0) always

**URL:** <https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757>\
**Category:** New to Julia\
**Created:** [March 10, 2025, 1:12pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757 "2025-03-10T13:12:43Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 10, 2025, 1:12pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/1 "2025-03-10T13:12:43Z")

</div>

Hello, I am currently trying to solve an ode where I have u(t) as a vector of length 1024, and I have du/dt = -u + A_M_f(u), where M is a matrix of 1024x1024. I specify some random u0, and solve the problem. I get a solution, but when I plot u(0) against u(100), as I am solving up to t = 100, I get the initial solution again. This happens even for simple problems, such as du = -u. Could someone help to explain what is going wrong? I have added the code below. I tried to make it self-contained, but if I made a mistake please let me know.

```julia
using DifferentialEquations
using Plots
pl = Plots

## Work to set up matrix 

function W(x)
    1/sqrt(π) .* exp.(-x.^2) .- 1/(sigma*sqrt(π)) .* exp.(-(x./sigma).^2)
end

n = 2^10; L = 10*π; h = (2*L)/n; 
x = -L .+ (0:n-1)*h
M = zeros(n,n);
y = h .* W(x)

start = n/2 + 1;

for i = 1:n
    M[i, :] .= circshift(y, start + i - 2)
end

#### Directly specifying initial conditions and parameters
σ = 1.5
A = 1.7
μ = 10.0
θ = 0.5

p = [σ, A, μ, θ]

u0 = 0.3*rand(n,1)

## Set up firing rate

function f(x)
    1/(1 .+ exp.(-mu .* x + theta)) - 1/(1 .+ exp(theta))
end

# Setting up the problem
@. rhs(u,p,t) = -u + A*M*f(u)

# Setting up trial problem

# @. rhs(du, u, p, t) = -u

# Time steps
tspan = (0.0,100.0)

## Solving problem

prob = ODEProblem(rhs, u0, tspan, p)

sol = solve(prob, reltol = 1e-6, saveat = 0.001)

final = sol.u[end]

p1 = pl.plot(final)
p2 = pl.plot(u0)
pl.plot(p1,p2, layout = (1,2))

```

---

<div class="post-metadata">

**Author:** ![Rohan\_Rajagopal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rohan_rajagopal/32/34183_2.png) [@Rohan\_Rajagopal](https://discourse.julialang.org/u/Rohan_Rajagopal)\
**Post date:** [March 10, 2025, 1:25pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/2 "2025-03-10T13:25:35Z")

</div>

Hi, I think you haven’t defined some of the functions(?)

```julia
W(x), WHat()

```

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [March 10, 2025, 1:31pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/3 "2025-03-10T13:31:34Z")

</div>

As @Rohan_Rajagopal said, your example code doesn’t quite work. You can open a fresh REPL and execute the code yourself and fix all the errors about undefined variables and then update it here.

But also without running the code, I’m a bit confused about the equations you’re trying to solve. `M` is a matrix, but then `f(u)` should be a vector and `M * f(u)` the regular matrix-vector product I guess?

Then your `rhs(u, p, t)` should probably look like (assuming `A` is just a number; note the missing dot between `M` and `f(u)`):

```julia
rhs(u, p, t) = - u .+ A .* M * f.(u)

```

(EDIT: I accidentally put a `*` instead of `+`)

and `f(x)` should have no dots, i.e. just map scalars to scalars.

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 10, 2025, 1:31pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/4 "2025-03-10T13:31:55Z")

</div>

> [@matyaspayn](#):
>
> `final = sol.u[end]`

Sorry, yes. Let me edit the original post

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 10, 2025, 1:36pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/5 "2025-03-10T13:36:04Z")

</div>

The code should work now. I thought since I put @. in front of my rhs I do not need to add the dots for the multiplication. The equation is -u (n,1 vector) + a scalar \* nxn matrix \* nx1 vector. The product of the matrix and f gives a vector, so Im just adding these vectors together with the scalar multiplication as well. I will try what you wrote, thanks

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 10, 2025, 1:38pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/6 "2025-03-10T13:38:52Z")

</div>

Unfortunately the issue persists. I should also note that since u is a vector, shouldn’t f have the dots, as its input is a vector? Or is this remedied by the f.(u)?

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [March 10, 2025, 1:41pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/7 "2025-03-10T13:41:30Z")

</div>

Yes, that’s almost correct, but the issue is that `@.` adds the dot to _every_ function call in the expression. But the matrix-vector product `M * x` has no dot, it’s just `*`.

Writing `M .* x` for a NxN matrix and a Nx1 vector will result in something different:

```julia-repl
julia> M = rand(3,3)
3×3 Matrix{Float64}:
 0.160868 0.997983 0.681036
 0.834087 0.619388 0.646875
 0.403048 0.658366 0.372333

julia> x = [1.0, 2.0, 3.0]
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> M * x
3-element Vector{Float64}:
 4.199942527126877
 4.013490057872009
 2.8367784827159843

julia> M .* x
3×3 Matrix{Float64}:
 0.160868 0.997983 0.681036
 1.66817 1.23878 1.29375
 1.20914 1.9751 1.117

```

This may or may not explain what you are seeing.

> [@matyaspayn](#):
>
> I should also note that since u is a vector, shouldn’t f have the dots, as its input is a vector? Or is this remedied by the f.(u)?

The latter. You can of course write the function `f(u)` such that it takes a vector and operates on it (with the dots as in your current example), but it’s common practice and in my opinion also easier to write and read the code if you write functions mapping scalars and then broadcast when you need it. The base functions also do it like that, e.g.

```julia-repl
julia> sin([1, 2, 3])
ERROR: MethodError: no method matching sin(::Vector{Int64})

julia> sin.([1, 2, 3])
3-element Vector{Float64}:
 0.8414709848078965
 0.9092974268256817
 0.1411200080598672

```

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 10, 2025, 1:49pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/8 "2025-03-10T13:49:55Z")

</div>

Okay thanks for help. I changed the code and checked that M\*f.(u) was an Nx1 matrix, which it was. But unfortunately the issue still persists. I will check to see if I have made mistakes like this elsewhere. Out of curiosity, in the example you showed, it says M \* x is a 3-element Vector. Mine says Matrix, and when I check the size it says Nx1. Is this a problem that it sees it as a matrix instead of a vector?

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [March 10, 2025, 3:39pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/9 "2025-03-10T15:39:17Z")

</div>

Yes, a `Vector{Float64}` and a `3×1 Matrix{Float64}` are different objects in Julia, but usually both represent what we would just call “vector” in maths. They might behave similarly, but they also have some important differences in code.

E.g. `*` again computes the “normal” matrix-vector product and spits out another vector, whereas `.*` broadcasts the vector over all columns of the matrix:

```julia-repl
julia> M = rand(3,3)
3×3 Matrix{Float64}:
 0.81235 0.683495 0.260397
 0.760464 0.447985 0.302931
 0.354402 0.05534 0.79264

julia> x = [1.0, 2.0, 3.0]; y = reshape(x, (3, 1))
3×1 Matrix{Float64}:
 1.0
 2.0
 3.0

julia> M * y
3×1 Matrix{Float64}:
 2.960533034374361
 2.5652254115093327
 2.843002725837259

julia> M .* y
3×3 Matrix{Float64}:
 0.81235 0.683495 0.260397
 1.52093 0.895969 0.605861
 1.06321 0.16602 2.37792

```

Note how the result of `M .* y` is again a matrix with each column multiplied (element-wise) with `[1,2,3]`.

* * *

Regarding your original question,

> [@matyaspayn](#):
>
> I get a solution, but when I plot u(0) against u(100), as I am solving up to t = 100, I get the initial solution again. This happens even for simple problems, such as du = -u.

If you ran the commented-out code from your post for `du = -u` then you will not see any change in your state because you are not updating `du` in your function. There are two different ways of specifying the ODE right-hand-side function and you mixed them up I believe. Here is an example of both versions side-by-side:

> **[Code Optimization for Differential Equations · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/dev/tutorials/faster_ode_example/)**
>
> Documentation for DifferentialEquations.jl.

What the actual code does or should do is still not completely clear to me, as it doesn’t run in a REPL (the parameters are sometimes called `sigma` and sometimes `σ` and they are defined after they are used. I tried to fix the issues below (also removing all the dots I thought are unnecessary and using a plain `Vector` instead of a `nx1 Matrix`) and I don’t get a constant solution for `u(t)`. But I have no idea if I’m solving the equations you intended to solve.

> **new code**
>
> ```julia
> using DifferentialEquations
> using Plots
> pl = Plots
> 
> ## Work to set up matrix 
> 
> σ = 1.5
> A = 1.7
> μ = 10.0
> θ = 0.5
> 
> n = 2^10
> L = 10*π
> h = (2*L)/n; 
> 
> p = [σ, A, μ, θ]
> 
> function W(x)
> 1/sqrt(π) * exp(-x^2) - 1/(σ*sqrt(π)) * exp(-(x/σ)^2)
> end
> 
> x = -L .+ (0:n-1) .* h
> M = zeros(n,n);
> y = h .* W.(x)
> 
> start = n/2 + 1;
> 
> for i = 1:n
> M[i, :] .= circshift(y, start + i - 2)
> end
> 
> #### Directly specifying initial conditions and parameters
> 
> u0 = 0.3*rand(n)
> 
> ## Set up firing rate
> 
> function f(x)
> 1/(1 + exp.(-μ * x + θ)) - 1/(1 + exp(θ))
> end
> 
> # Setting up the problem
> rhs(u,p,t) = -u + A .* M * f.(u)
> 
> # Setting up trial problem
> 
> #### This function doesn't work as expected, since you would need to do `du .= -u`
> # @. rhs(du, u, p, t) = -u
> 
> # Time steps
> tspan = (0.0,1.0)
> 
> ## Solving problem
> 
> prob = ODEProblem(rhs, u0, tspan, p)
> 
> sol = solve(prob, reltol = 1e-6, saveat = 0.1)
> 
> final = sol.u[end]
> 
> p1 = pl.plot(final)
> p2 = pl.plot(u0)
> pl.plot(p1,p2, layout = (1,2))
> 
> ```

* * *

PS: The manual has a section about array arithmetic and broadcasting that might be helpful:

> **[Single- and multi-dimensional Arrays · The Julia Language](https://docs.julialang.org/en/v1/manual/arrays/#man-array-and-vectorized-operators-and-functions)**
>
> Documentation for The Julia Language.

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 10, 2025, 5:47pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/10 "2025-03-10T17:47:09Z")

</div>

Thank you for the links. I will take a look. I am using JupyterNotebook, so maybe this is the issue with the REPL. I tried the code you wrote but unfortunately I still get u(end) = u(0). I also noticed that when I do M \* f.(0) I still get a matrix instead of a vector. Not sure how to remedy this, as in the code you sent earlier M \* x was a 3-element vector,

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 11, 2025, 4:30am UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/11 "2025-03-11T04:30:44Z")

</div>

> [@matyaspayn](#):
>
> I also noticed that when I do M \* f.(0) I still get a matrix instead of a vector.

Yes because `f.(0)` is a scalar. So `M*f.(0)` just multiplies the matrix `M` by a scalar which yields another matrix.

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 11, 2025, 4:43pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/12 "2025-03-11T16:43:10Z")

</div>

Dear all, thank you for the help. I logged in today and restarted my kernel in JupyterNotebook and now the code works. I have had this happen a couple times now, where I change nothing and code goes from either working to not or vice versa, but hopefully this does not happen again.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 11, 2025, 7:42pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/13 "2025-03-11T19:42:04Z")

</div>

> [@matyaspayn](#):
>
> I have had this happen a couple times now, where I change nothing and code goes from either working to not or vice versa, but hopefully this does not happen again.

That is an inherent weakness of Jupyter notebooks… Have you tried Pluto.jl notebooks? They’ll always keep your notebook in a consistent state. I cannot recommend them enough! 🙂

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 11, 2025, 8:09pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/14 "2025-03-11T20:09:12Z")

</div>

Just tried it and loving it. Thanks for the recommendation!

One thing: I noticed that my output always shows, even if I put a ; at the end. Any quick fix to this?

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 11, 2025, 8:14pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/15 "2025-03-11T20:14:43Z")

</div>

Strange. I think usually a ; should suppress the output. If it’s important to you, you cann add `; nothing` to explicitly return `nothing`

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 11, 2025, 8:17pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/16 "2025-03-11T20:17:08Z")

</div>

Also another recommendation: There is PlutoUI.jl which adds a couple of nice features on top. See the examples here:

> **[Interactivity](https://featured.plutojl.org/basic/plutoui.jl)**
>
> Sliders, buttons, dropdowns, and more using PlutoUI.jl!

Especially the `TableOfContents` can be found in virtually all of my notebooks 🙂

---

<div class="post-metadata">

**Author:** ![matyaspayn](https://avatars.discourse-cdn.com/v4/letter/m/57b2e6/32.png) [@matyaspayn](https://discourse.julialang.org/u/matyaspayn)\
**Post date:** [March 12, 2025, 9:16am UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/17 "2025-03-12T09:16:27Z")

</div>

Alright this works. Thanks! Another issue is if i define lets say A = 1 in one cell and then update the value of A later in a different cell, writing A = 2, this causes an error. Is this because Pluto.jl needs a different way of updating variables?

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 12, 2025, 12:19pm UTC](https://discourse.julialang.org/t/ode-returning-f-t-f-0-always/126757/18 "2025-03-12T12:19:28Z")

</div>

Yes this is what sets Pluto.jl aparat from usual notebook systems. You will need to adapt your habits a bit.  
Ask yourself: Why do you want to redefine `A`?

- is it conceptually the same `A`? Then perhaps you should change `A` in the cell you currently define it and Pluto will run all cells that depend on `A` automatically for you.
- is it conceptually a different `A` because you want to repeat some computation with a different value? Then I recommend choosing a more descriptive variable name to distinguish your two `A`s

With Pluto I am generally using longer and more descriptive global variable names. Which maybe makes it slightly more annoying to write code, but much much easier reading the code later 🙂