# Solving and ODE system starting from different initial points

Hello,

I have a simple system of which i would like to build the phase portrait.
I thought of creating a mesh of interesting points and solve the problem for each of them.
This is what i naively tried until now but it does not work.

xs = -0.5:0.1:1.5
ys = -0.5:0.1:3
u0s = [[x,y] for x in xs, y in ys]
tspan = (0.0, 200.0)
λ = 0
prob = ODEProblem(f, u0s, tspan, λ)
sol = solve(prob, progress=true)


I could not find anything specific in the documentation on how to approach this problem, i know about remake but i am not sure if it is useful or a good idea to use it with so many new initialisations.

Could a for loop work where every time i store the solution for every initial point?
something on the likes of

for u0 in u0s
sol = solve(remake(prob; u0=u0))
end


and then store sol somewhere for later access and/or plotting.
(in this case, how do i store the solutions correctly?)

Any smarter ideas, this feels really dumb.

Thanks!

Two possibilities:

• Shouldn’t for u0 in eachindex(u0s) instead be for u0 in u0s?
• The remake documentation suggests you should use keyword assignment; this may be unnecessary, though?? Thus sol = solve(remake(prob; u0 = u0))

You are right about both points, i fixed the code and now at least runs.
But now if i try to do

plot!(sol)


with the new computed solution it does not get updated.

Anyhow, i was mainly wondering if the ODEProblem interface accepted some kind of “multiple initial point” kind of thing, the current approach feels rather hacky.

The following works – I’ve used a boring pendulum model, since you don’t give your model:

function pendulum(du,u,p,t)
g = 9.81
L = 2
alpha = u
omega = u
du = omega
du = -g/L*sin(alpha)
end
#
as = -2pi:pi/6:2pi
ws = -0.5:0.1:0.5
us = [[a,w] for a in as, w in ws]
tspan = (0,5.)
#
prob = ODEProblem(pendulum,us,tspan)
sol = solve(prob)
plot(sol,vars=(1,2),label="")
#
for u in us
sol = solve(remake(prob,u0 = u))
plot!(sol,vars=(1,2),label="")
end
#
plot!()


The result is as follows (not very interesting…):

Thanks a bunch!
I had tried before without the final plot!() and it did not work since i forgot the display inside the for loop!

1 Like

I had trouble finding the right terminology in the documentations, thanks for pointing me to it!