Hi,

I have the sort of thankless job of writing up a Julia-implementation for what I guess is a matrix-version of a Riccati equation using the DifferentialEquations.jl package. More spesifically this is a boundary value problem.

I think I ran into some problems regarding complex numbers, and type conversion - and could use some pointers on type conversion etc.

More spesifically, I think I need help type hinting the ComplexF64 type everywhere. I have a hunch it is the formatting of the initial guess that is off. The documentations says I could supply it as a function of the free variable, here x, but not how to do so. The way I have done it now, I have 16 coupled equations, but I mathematically speaking I only need 4 coupled equations for 2x2 matrices.

```
#Functions
equation!(dg, g, p, x) #This encodes the equation - with the independent variable being x and not t.
bc!(res, sol, p, x) #Again, sorry, x is the independent variable.
#Calls
bvp1 = BVProblem(equation!, bc!, g_0, xspan)
sol1 = solve(bvp1, Shooting(Vern7()), dt = 0.03) #dt should be read dx
```

Can you share an example code which demonstrates your issue?

Sure! I was a bit hesitant to do so - as I am probably breaking every best-practice imaginable. Here is a simpler version with scalar values. Chances are I created some numerical instabilities now.

EDIT: Having the initial guess as zeros produces a valid solution, and as ones a singularity - now it does not do that

```
function simple_eqn!(dg, g, p, x)
#Here, I think I want dg and g to be Vector{ComplexF64}, but they look real
γ1, γ2, dγ1, dγ2 = g[1:4]
N = 1/(1 - γ1*γ2) #Readability
dg[1] = dγ1
dg[2] = dγ2
dg[3] = - 2 * dγ1 * γ2 * N * dγ1 - 2im * E * γ1
dg[4] = - 2 * dγ2 * γ1 * N * dγ2 + 2im * (-E) * γ2
end
function simple_bc!(res, sol, p, x)
res[1] = sol[1][3]
res[2] = sol[1][4]
res[2] = sol[end][3]
res[2] = sol[end][4]
#Chances are I did this wrong, but I want the derivatives, i.e. dγ1, dγ2 to be zero at the edges
end
E = 0.2 #Physical constant
g_0 = 0.5 * ones(ComplexF64, 4)
xspan = (0.0, 1)
bvp1 = BVProblem(simple_eqn!, simple_bc!, g_0, xspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt = 0.03) #dt should be read dx
```