For the special case of unitary matrices, my quantum mechanics books tells me to use Time Evolution Operators.

`β(t) = exp(-im*H*t/h)*β(0)`

As far I know, I just need to use Matrix Exponential.

First, decompose `H = ψ*λ*inv(ψ)`

Second, we create a time evolution operator `U = exp(-im*H*t) = ψ*exp(-im*λ*t)*inv(ψ)`

If one sets a small enough time interval `δt`

: `δU = exp(-im*H*δt) = ψ*exp(-im*λ*δt)*inv(ψ)`

Then, I can evolve iteratively the solution `βₜ₊₁ = δU*βₜ`

of my problem.

```
using LinearAlgebra
using Random
Random.seed!(1256)
N = 2
H = [1 1;-1 1]./sqrt(2)
β₀= rand(ComplexF64, N);
maxTimeSteps, t_max = 1000, 10
t = range(1, stop=t_max, length=maxTimeSteps);
δₜ = t[2]-t[1];
@time begin
λ, ψ = eigen(H);
δU = ψ*exp(-im*Diagonal(λ*δₜ))*inv(ψ)
βₜ = []
push!(βₜ, deepcopy(β₀));
for idx_t = 2:maxTimeSteps
push!(βₜ, δU*βₜ[idx_t-1] )
end
end
```

I would like to solve this problem with `DifferentialEquations`

, but I don’t know how to properly write down. It should be simple, but I’m taking time.

My bet was to use `du/dt = a*u`

as input, because its solution is `u(t)=exp(a*t)`

```
using DifferentialEquations
function ff(du, u, p,t)
du[:] .= -im*H*u
end
tspan = (0.0,t_max);
prob = ODEProblem(ff,β₀,tspan);
@time sol = solve(prob, Tsit5(), reltol=δₜ, abstol=δₜ);
```

I can check the difference of the last results are far away of each other.

```
julia> abs.(βₜ[end] .- sol.u[end])
2-element Vector{Float64}:
424.06945981919637
424.06825305239903
```

Is there any one with physics background that already knows where my silly mistake is ?