Have you tried having a look a the simulation output provided by the double_pendulum.jl code that you are referencing here? In particular, from the first simulation you get

```
julia> sol1
retcode: Success
Interpolation: automatic order switching interpolation
t: 134-element Array{Float64,1}:
0.0
0.005
0.05499999999999999
⋮
24.71989277613792
24.86868170266417
25.0
u: 134-element Array{Array{Float64,1},1}:
[1.0, 1.0, 1.0, 1.0]
[0.9173873470403654, 0.9174516957353678, 1.0023967042847444, 1.0000000563149365]
[0.08644609908203948, 0.09196004700351201, 1.0149015079311146, 1.0000567812684389]
⋮
[-6.054231992074114, -12.056668232377719, 0.027646098435418, 0.10612941900049216]
[-6.1495227859460195, -11.69036787508306, -6.080301318376474e-5, -0.32640770884302933]
[-6.256313663905104, -10.239601242478551, -0.1324936832683153, -0.629101072316051]
```

Obviously, the simulation output comes as some more complicated data type than just a an array (in fact, a matrix) that you are used to from Matlab. But getting the array of values is straigtforward using the dot. Note that the vector of (state) variables is denoted `u`

here:

```
julia> sol1.u
134-element Array{Array{Float64,1},1}:
[1.0, 1.0, 1.0, 1.0]
[0.9173873470403654, 0.9174516957353678, 1.0023967042847444, 1.0000000563149365]
[0.08644609908203948, 0.09196004700351201, 1.0149015079311146, 1.0000567812684389]
[-1.9803904064051525, -1.9817102851965904, 0.9544687664602022, 1.0006037605753557]
[-4.5565663813837265, -4.95682194749687, 0.672002327426011, 0.9832650510821658]
[-5.9032809983798655, -7.790178097454042, 0.3242060850308398, 0.8759059209509462]
[-5.943613484795805, -10.775076914267512, 0.022043544720056207, 0.5667541455845946]
[-5.7130796782664195, -12.013371802990793, -0.0426856254161595, 0.14784486944733183]
[-5.648958300778203, -11.889774751071544, -0.007217148908584846, -0.25005242797974475]
⋮
[5.587060546627331, 7.806700961000203, 0.37552602998220136, 0.8960073447463703]
[3.8453217414924703, 4.448148406573896, 0.7144888073156342, 1.0376447446638644]
[1.286466475379462, 1.376402400305727, 0.9166785732656021, 1.0694682774712758]
[-2.491704854444083, -2.7950277723838823, 0.8532718060615504, 1.0576089418582855]
[-5.274180306238754, -6.898084330762543, 0.47881885731568813, 0.9453633766685099]
[-6.015148096215678, -10.687330502172246, 0.11381677838517706, 0.5827117250218716]
[-6.054231992074114, -12.056668232377719, 0.027646098435418, 0.10612941900049216]
[-6.1495227859460195, -11.69036787508306, -6.080301318376474e-5, -0.32640770884302933]
[-6.256313663905104, -10.239601242478551, -0.1324936832683153, -0.629101072316051]
```

Note that it is actually an array of arrays, the outer dimension corresponding to time, therefore, the initial state vector is obtained with

```
julia> sol1.u[1]
4-element Array{Float64,1}:
1.0
1.0
1.0
1.0
```

Now, obtaining a function of this sequence of state vectors applied *time-wise* is easy. Perhaps you may want to get the following function

```
julia> f(u) = cos(u[1]) + u[2]*u[3] + 5*u[4]
f (generic function with 1 method)
```

Applying this to the simulation outcomes is straightforward using broadcasting:

```
julia> f.(sol1.u)
134-element Array{Float64,1}:
6.54030230586814
6.527547555276698
6.089880158973252
2.7133012143705875
1.4301365857910806
2.7826065872357417
3.5391472089129405
2.0938666057229924
-0.358921886704315
⋮
8.178989078710751
7.603941442605557
6.889574248219444
2.106974617987561
1.9565870812663815
2.6614536749009683
1.1712317404892454
-0.6402472780756148
-0.7891838986468183
```

Is this what you wanted?