I want to write gradient ascent pulse engineering (GRAPE) and I try the simple Hamiltonian H=\Omega*|0\rangle\langle 1| + \Omega*|1\rangle\langle 0|. The initial state is |0\rangle\langle 0| and the final state is |1\rangle\langle 1|. However, I cannot get the right result. I have two questions: (1) I cannot show my \Omega pulse in figure; (2) The GRAPE dose not achieve the desired outcome. I cannot find what my question is. Here is my code:

```
# "gradient ascent pulse engineering (GRAPE) coed demo"
using LinearAlgebra, PyPlot
# Define variables
dim = 2
Natom = 1
# Define basis
s0 = [1; 0]
s1 = [0; 1]
tr01 = s0 * s1'
tr10 = s1 * s0'
# Start state
global Ļ0 = s0 * s0'
# Target state
global Ļf = s1 * s1'
# Time list
t_list = range(0, 1, length=101)
dt = t_list[2] - t_list[1]
# Initial Pulses
global Ī© = rand(length(t_list)) .* (2Ļ) .- Ļ
# Initialize
fidelity = []
Ļ0_list = Vector{Any}(undef, length(t_list))
Ļf_list = Vector{Any}(undef, length(t_list))
grad_Ī© = zeros(length(t_list))
function PWC_pulse(pwc_pulse)
pwc_pulse = [0; pwc_pulse; 0]
time_steps = 0:(length(pwc_pulse)-1)
time_steps_stair = repeat(time_steps, 1, 2)[1:end-1]
pwc_pulse_stair = repeat(pwc_pulse, 1, 2)[1:end-1]
return pwc_pulse_stair
end
for iter in 1:1501
# Bidirectional Evolution
for i in 1:length(t_list)
H = Ī©[i] * (s0 * s1' + s1 * s0')
U_indt = exp(-im * H * dt)
global Ļ0 = U_indt * Ļ0 * U_indt'
Ļ0_list[i] = Ļ0
end
for i in length(t_list):-1:1
H = Ī©[i] * (s0 * s1' + s1 * s0')
U_indt = exp(im * H * dt)
global Ļf = U_indt' * Ļf * U_indt
Ļf_list[i] = Ļf
end
# Gradient Calculation
for i in 1:length(t_list)
grad_Ī©[i] = -real(tr(Ļf_list[i]' * 1im * dt * ((tr01 + tr10) * Ļ0_list[i] - Ļ0_list[i] * (tr01 + tr10))))
end
# Update Pulses
lr = 0.01 - iter / 100
global Ī© = Ī© + lr * grad_Ī©
Ī© = max.(-Ļ, min.(Ļ, Ī©))
# Fidelity
global fidelity = abs(tr(Ļf * (s1 * s1')))
global infid = 1 - fidelity
# Plot Pulses
if mod(iter, 100) == 0
pwc_pulse_stair = PWC_pulse(Ī©)
figure()
plot(pwc_pulse_stair)
xlabel("time")
ylabel("pulse strength")
title("$(iter)-th fidelity: $(round(fidelity, digits=4))")
if iter == 300
savefig("pulses_300.png")
end
sleep(0.01)
end
end
figure()
plot(1:length(infid), log.(infid))
xlabel("iteration")
ylabel("infidelity")
savefig("infidelity.png")
gcf()
```