I wouldn’t use the phrase “evaluate a *symbolics* variable”. At this stage, the solutions are available and they are *numeric*. But anyway, maybe your phrase is ok.

It is important to understand the *structure* of the solution which I have called `sol`

. The solution contains two fields `t`

and `u`

– field `t`

holds the time instances where the solution has been found, while field `u`

holds the computed values of the unknowns in these time instances.

In the case of Backstress:

```
> sol.t' # I use transpose `'` symbol to save vertical space
1×32 adjoint(::Vector{Float64}) with eltype Float64:
0.0 0.762987 8.39286 84.6916 … 44162.9 45314.9 47666.0 50000.0
> sol.u'
1×32 adjoint(::Vector{Vector{Float64}}) with eltype LinearAlgebra.Adjoint{Float64, Vector{Float64}}:
[1.0;;] [1.00003;;] [1.00033;;] [1.00329;;] … [70.7544;;] [70.751;;]
```

Observe that we have found the solution in 32 time instances. Also, observe that `sol.u`

is a *vector of vectors*. Each element in the vector holds the solution of all unknowns in the corresponding time instance. So as an example, vector element `sol.u[3]`

holds the solution of unknowns corresponding to time `sol.t[3]`

. Two specific points:

- Why is there only one unknown with solution in
`sol.u`

?
- The shape of
`sol.u`

is not in a form that can be plotted.

`sol.u`

and unknowns.

Field `u`

of the solution only holds the solution to variables formally defined as `unknowns`

. To find out which these are, and their order of storage in `sol.u`

, you need to query:

```
> unknowns(sys)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
σ_T(t)
```

What about other “unknowns”? These have been stripped off the model and made into `observed`

variables; the observed variables are found by:

```
> observed(sys) |> println
Equation
[ϵ°(t) ~ A°*Exp*
sinh((-σ_T(t) + σ - σ_d) / σ_ref)*(σ_T(t)^2)]
```

I’ll get back to how to compute this observed variable – which is what you ask about, I guess.

- The shape of
`sol.u`

and how to make it plottable…

Perhaps the most efficient way to make `sol.u`

plottable is:

```
> reduce(vcat,sol.u)' # Again: transpose to save space
1×32 adjoint(::Vector{Float64}) with eltype Float64:
1.0 1.00003 1.00033 1.00329 … 70.7375 70.7376 70.7544 70.751
```

Now you have the solution in a vector of elements, and not in a vector of vectors. The same should work if `sol.u`

has multiple unknowns.

Plotting the result can be done by:

```
plot(sol.t, reduce(vcat,sol.u))
```

which produces:

Notice how jagged the solution is compared to doing:

```
plot(sol)
```

which produces

The reason for the jaggedness in the first method, `plot(sol.t, reduce(vcat,sol.u))`

, is that this draws straight lines between the datapoints of `sol.t[j]`

and `sol.u[j]`

.

One of the fancy things with Julia solvers is that `sol`

not only holds the solution in datapoints, i.e., `sol.t`

, etc. The solution also *includes* interpolation functions. Interpolation in a certain time point `t0`

is given by syntax `sol(t0)`

, e.g.:

```
> sol(4e4)
1-element Vector{Float64}:
70.74021605315112
```

This is the solution ` σ_T(t=4e4)`

.

Final question: how can you *specify* which variable you want to know the value of? And how can you find the value of *observed* variables? Simple:

```
> sol(4e4, idxs=σ_T)
70.74021605315112
> sol(4e4, idxs=ϵ°)
3.6413814395589957e-7
> sol(4e4, idxs=[σ_T, ϵ°])
2-element Vector{Float64}:
70.74021605315112
3.6413814395589957e-7
```

What about plotting observed variables (and unknowns) with interpolation?

```
plot(sol, idxs=σ_T)
```

produces.

It is also possible to scale time variable and results:

```
plot(sol, idxs=(t/60,[ϵ°*1e8,σ_T]), label=["epsilon*1e8" "sigma"], xlabel="time [min]", ylabel="")
```

produces:

If you want a parametric plot of the two dependent variables with time as parameter:

```
plot(sol, idxs=(σ_T, ϵ°*1e8))
```

produces: