I am using `ODEProblem`

to solve a system of ODEs: https://mathb.in/75542 such that u_1 = u_N = 0.

To do this, I wrote the code shown at the end of the post to generate this solution: https://i.imgur.com/wQC8aTj.png I expected the solution to decay near the boundary, but instead we have quite the opposite.

**Question: ** Does the following code indeed solve the system (1) shown in https://mathb.in/75542 ?

```
using DifferentialEquations, GLMakie, LinearAlgebra, SpecialFunctions
GLMakie.activate!()
#=
DifferentialEquations - Used for the ODE Solver
GLMakie - 2D/3D plotting backend,
LinearAlgebra,
SpecialFunctions - Used for importing the Riemann zeta function
=#
const N = 31
centerGrid = (N-1)/2;
x = -centerGrid:centerGrid;
#f(n,α) = sum(g(n,m,α) for m in 1:N if m!=n) # Diagonal elements for the coefficient matrix M
#g(n,m,α) = 1/abs(n-m)^(1+α) # Off-diagonal elements for the coefficient matrix M
function main()
p = [1] # ODEProblem() requries a parameter list: `p`. Here, p = [h]
α = 0.2
h = 1
tspan = (0.0,100)
RHS(u,p,t) = 1*im*(p[1]*coefficientMatrix(α)*u - @.(((abs(u))^2).*u) )
initial_values = Complex.(sech.(x))
# Equivalent to Ode45 from Matlab
prob = ODEProblem(RHS, initial_values, tspan, [h])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
plot_main(sol) # Surface and contour plot of |u|²
return sol
end
function diagonal_coefficients(n,α)
if n == 1 || n == N
diag_value = 0
else
diag_value = 1/abs(n-1)^(1+α) + sum(1/abs(n-m)^(1+α) for m in 2:N-1 if m!=n) + 1/abs(n-N)^(1+α) # Entry for M[n,n]. row n, col n
end
return diag_value
end
function off_diagonal_coefficients(n,m,α)
if n == 1 || n == N
off_diag_value = 0
else
off_diag_value = -1/abs(n-m)^(1+α) # Entry for M[n,m]. row n, col m
end
return off_diag_value
end
function coefficientMatrix(α)
# M is the coefficient matrix corresponding to the 1-dimensional nonlinear discrete operator of the fDNLS
M = Matrix{Float64}(undef, length(x), length(x))
for i in 1:size(M,1) # Iterate through the rows
for j in 1:size(M,2) # Iterate through the columns
if i==j
M[i,i] = diagonal_coefficients(i, α)
elseif i!=j
M[i,j] = off_diagonal_coefficients(i,j,α)
end
end
end
return M
end
function plot_main(sol)
fig = Figure(backgroundcolor=:snow2)
y = sol.t
z= [abs(sol[i,j])^2 for i=1:size(sol)[1], j=1:size(sol)[2]] # sol[i,j] is the ith component at timestep j. # sol[i,j] is the ith component at timestep j.
p1 = surface(fig[1,1], x, y, z,axis=(type=Axis3, xlabel="Space", ylabel="Time", zlabel="|u|²", xticks=-centerGrid:2:centerGrid))
p2 = contourf(fig[1,2], x, y, z, axis=(xlabel="Space", ylabel="Time"))
display(fig)
end
main()
```