Unexpected behavior when returning DifferentialEquations solve() solution from external .jl file

I am new to Julia. I am trying to use the DifferentialEquations.jl to solve, well, some ODEs, via Python. I can, for example, use the following python script to solve the Lorenz attractor:

from diffeqpy import de
from julia import Main

jul_f = Main.eval("""
function f(dy,y,p,t)
    x, y, z = y
    sigma, rho, beta = p
    dy[1] = sigma * (y - x)
    dy[2] = x * (rho - z) - y
    dy[3] = x * y - beta * z
end""")
u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,2.66]
prob = de.ODEProblem(jul_f, u0, tspan, p)
sol = de.solve(prob)

after which sol is an object with attributes I can interrogate, e.g. sol.t contains the t values and sol.u contains the solution.

In comparison, if I put everything in a separate file, MySol.jl like this:

using DifferentialEquations
using BenchmarkTools

function f(dy,y,p,t)
    x, y, z = y
    sigma, rho, beta = p
    dy[1] = sigma * (y - x)
    dy[2] = x * (rho - z) - y
    dy[3] = x * y - beta * z
end

function solver(u0, t_start, t_end, p)
    problem = ODEProblem(f, u0, (t_start, t_end), p)
    sol = solve(problem)
    return sol
end

and then in python do the following:

from julia import Main

Main.include("MySol.jl")

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,2.66]

sol = Main.solver(u0, t_span[0], t_span[1], p)

then sol is simply an array which contains the (correct) y-values for the solution and nothing of the remaining values which would be in the struct. Is there a way to fix this so that I can recover the whole solution object and not just the y-values?

@c42f maybe you might know why PythonCall would convert here? I don’t know why it would do a conversion if you use include like that.