Jump

    # State Space
    n = length(x)

    # Initialize Model
    m = Model(solver = IpoptSolver(max_iter = iterations, constr_viol_tol = maxdist))

    # Parameters
    @NLparameter(m, ρ == byp.ρ)
    @NLparameter(m, δ == byp.δ)
    @NLparameter(m, B1 == byp.B1)
    @NLparameter(m, δ1 == byp.δ1)
    @NLparameter(m, B2 == byp.B2)
    @NLparameter(m, δ2 == byp.δ2)
    @NLparameter(m, ω == byp.ω)
    @NLparameter(m, μ == byp.μ)
    @NLparameter(m, σ == byp.σ)
    @NLparameter(m, νA == byp.νA)
    
    @NLparameter(m, γA == byp.γA)
    @NLparameter(m, ψA == byp.ψA)
    @NLparameter(m, γB == byp.γB)
    @NLparameter(m, ψB == byp.ψB)

    # Variables to minimize
    @variable(m, gA[i=1:n], start = initial[0 * n + i])
    @variable(m, gB[i=1:n], start = initial[1 * n + i])
    @variable(m, ϕ1[i=1:n], start = initial[2 * n + i])
    @variable(m, ϕ2[i=1:n], start = initial[3 * n + i])

    # Expressions
    ## Constants
    @NLexpression(m, αA, 1 - 1/ψA)
    @NLexpression(m, αB, 1 - 1/ψB)
    @NLexpression(m, M1A, - 1 - (-(αA + γA -1) / αA) / γA)
    @NLexpression(m, M2A, (γA - 1) / γA)
    if byp.γB < Inf
        @NLexpression(m, M1B, - 1 - (-(αB + γB -1) / αB) / γB)
        @NLexpression(m, M2B, (γB - 1) / γB)
    end

    @NLexpression(m, invx, 1 / (x[2]-x[1]))
    ## First Derivative
    @NLexpression(m, gAp[i=1:n], (i == 1) * (- 0.5 * (3 * gA[1] - 4 * gA[2] + gA[3]) * invx) + (i == n) * (0.5 * (3 * gA[n] - 4 * gA[n - 1] + gA[n - 2]) * invx) + (i > 1) * (i < n) * (0.5 * (gA[min(i + 1, n)] - gA[max(i - 1, 1)]) * invx))
    @NLexpression(m, gBp[i=1:n], (i == 1) * (- 0.5 * (3 * gB[1] - 4 * gB[2] + gB[3]) * invx) + (i == n) * (0.5 * (3 * gB[n] - 4 * gB[n - 1] + gB[n - 2]) * invx) + (i > 1) * (i < n) * (0.5 * (gB[min(i + 1, n)] - gB[max(i - 1, 1)]) * invx))
    @NLexpression(m, ϕ1p[i=1:n], (i == 1) * (- 0.5 * (3 * ϕ1[1] - 4 * ϕ1[2] + ϕ1[3]) * invx) + (i == n) * (0.5 * (3 * ϕ1[n] - 4 * ϕ1[n - 1] + ϕ1[n - 2]) * invx) + (i > 1) * (i < n) * (0.5 * (ϕ1[min(i + 1, n)] - ϕ1[max(i - 1, 1)]) * invx))
    @NLexpression(m, ϕ2p[i=1:n], (i == 1) * (- 0.5 * (3 * ϕ2[1] - 4 * ϕ2[2] + ϕ2[3]) * invx) + (i == n) * (0.5 * (3 * ϕ2[n] - 4 * ϕ2[n - 1] + ϕ2[n - 2]) * invx) + (i > 1) * (i < n) * (0.5 * (ϕ2[min(i + 1, n)] - ϕ2[max(i - 1, 1)]) * invx))
    ## Second Derivative
    @NLexpression(m, gApp[i=1:n], (i == 1) * ((gA[3] + gA[1] - 2 * gA[2]) * invx^2) + (i == n) * ((gA[n] + gA[n - 2] - 2 * gA[n - 1]) * invx^2) + (i > 1) * (i < n) * ((gA[min(i + 1, n)] + gA[max(i - 1, 1)] - 2 * gA[i]) * invx^2))
    @NLexpression(m, gBpp[i=1:n], (i == 1) * ((gB[3] + gB[1] - 2 * gB[2]) * invx^2) + (i == n) * ((gB[n] + gB[n - 2] - 2 * gB[n - 1]) * invx^2) + (i > 1) * (i < n) * ((gB[min(i + 1, n)] + gB[max(i - 1, 1)] - 2 * gB[i]) * invx^2))
    @NLexpression(m, ϕ1pp[i=1:n], (i == 1) * ((ϕ1[3] + ϕ1[1] - 2 * ϕ1[2]) * invx^2) + (i == n) * ((ϕ1[n] + ϕ1[n - 2] - 2 * ϕ1[n - 1]) * invx^2) + (i > 1) * (i < n) * ((ϕ1[min(i + 1, n)] + ϕ1[max(i - 1, 1)] - 2 * ϕ1[i]) * invx^2))
    @NLexpression(m, ϕ2pp[i=1:n], (i == 1) * ((ϕ2[3] + ϕ2[1] - 2 * ϕ2[2]) * invx^2) + (i == n) * ((ϕ2[n] + ϕ2[n - 2] - 2 * ϕ2[n - 1]) * invx^2) + (i > 1) * (i < n) * ((ϕ2[min(i + 1, n)] + ϕ2[max(i - 1, 1)] - 2 * ϕ2[i]) * invx^2))
    ## Other obiects
    if byp.γB < Inf
        @NLexpression(m, Γ[i=1:n], 1 / (x[i] / γA + (1 - x[i]) / γB))
    else
        @NLexpression(m, Γ[i=1:n], 1 / (x[i] / γA))
    end
    @NLexpression(m, Θ[i=1:n], x[i] / (1 - αA) + (1 - x[i]) / (1 - αB))
    @NLexpression(m, ωA[i=1:n], x[i] * Γ[i] / γA)
    if byp.γB < Inf
        @NLexpression(m, ωB[i=1:n], (1 - x[i]) * Γ[i] / γB)
        @NLexpression(m, σX[i=1:n], x[i] * (Γ[i] - γA) / (Γ[i] / γB * x[i] * (1 - x[i]) * ((1 - γA - αA) / αA * gAp[i] / gA[i] - (1 - γB - αB) / αB * gBp[i] / gB[i]) + γA) * σ)
        @NLexpression(m, κ[i=1:n], Γ[i] * σ + ωA[i] * ((1 - γA - αA) / αA) * gAp[i] / gA[i] * σX[i] + ωB[i] * ((1 - γB - αB) / αB) * gBp[i] / gB[i] * σX[i])
        @NLexpression(m, nB[i=1:n], (2 - αB) / (2 * γB * (1 - αB)) * κ[i]^2 + (αB + γB - 1) / (2 * γB * αB) * (gBp[i] / gB[i] * σX[i])^2 + (αB + γB -1) /(γB * αB)  * (gBp[i] / gB[i] * σX[i]) * κ[i])
    else
        @NLexpression(m, σX[i=1:n], x[i] * (Γ[i] / γA - 1) * σ / (1 + x[i] * ((Γ[i] * x[i] / γA - 1) * (1 - ψA * γA) / γA * gAp[i] / gA[i] +  Γ[i] * (1 - x[i]) / γA * (1 - ψB) * gBp[i] / gB[i])))
        @NLexpression(m, κ[i=1:n], Γ[i] * (σ - x[i] * (1 - ψA * γA) / γA * gAp[i] / gA[i] * σX[i] - (1 - x[i]) * (1 - ψB) * gBp[i] / gB[i] * σX[i]))
        @NLexpression(m, nB[i=1:n], 1 / 2 * (gBp[i] / gB[i] * σX[i])^2 * (1 - ψB) * (γA - ψB))
    end
    @NLexpression(m, nA[i=1:n], (2 - αA) / (2 * γA * (1 - αA)) * κ[i]^2 + (αA + γA - 1) / (2 * γA * αA) * (gAp[i] / gA[i] * σX[i])^2 + (αA + γA -1) /(γA * αA)  * (gAp[i] / gA[i] * σX[i]) * κ[i])
    @NLexpression(m, r[i=1:n], ρ + 1 / Θ[i] * (μ - δ * (νA * gA[i] * (ϕ1[i] + ϕ2[i]) + (1 - νA) * gB[i] * (ϕ1[i] + ϕ2[i]) - 1)) - 1 / Θ[i] * (x[i] * nA[i] + (1 - x[i]) * nB[i]) )
    @NLexpression(m, μX[i=1:n], x[i] * ((r[i] - ρ)/(1 - αA) + nA[i] - δ - μ) + νA * δ * gA[i] * (ϕ1[i] + ϕ2[i]) - σ * σX[i])

   

    # Constraints
    @NLconstraint(m, gAPDE[i=1:n], (σX[i]^2 / 2 * M1A * ((M1A - 1) * gAp[i]^2 / gA[i] + gApp[i]) + M1A * gAp[i] * (μX[i] - M2A * σX[i] * κ[i]) + (κ[i]^2 / 2 * M2A * (M2A - 1) - M2A * (r[i] + δ) - M1A * gA[i] +  (-(ρ + δ) / αA * (1 - γA)) / γA) * gA[i]) / M1A == 0)
    if byp.γB < Inf
        @NLconstraint(m, gBPDE[i=1:n], (σX[i]^2 / 2 * M1B * ((M1B - 1) * gBp[i]^2 / gB[i] + gBpp[i]) + M1B * gBp[i] * (μX[i] - M2B * σX[i] * κ[i]) + (κ[i]^2 / 2 * M2B * (M2B - 1) - M2B * (r[i] + δ) - M1B * gB[i] +  (-(ρ + δ) / αB * (1 - γB)) / γB) * gB[i]) / M1B == 0)
    else
        @NLconstraint(m, gBPDE[i=1:n], σX[i]^2 / 2 * (gBpp[i] + gBp[i]^2 / gB[i] * (ψB / (1 - ψB) - γA / (1 - ψB))) +  gBp[i] * μX[i] + ((1 - ψB) * (r[i] + δ) + ψB * (ρ + δ) - gB[i]) * gB[i] == 0)
    end

    @NLconstraint(m, ϕ1PDE[i=1:n], ϕ1pp[i] * σX[i]^2 / 2 + ϕ1p[i] * (μX[i] + σX[i] * (σ - κ[i])) + ϕ1[i] * (μ - r[i] - δ - δ1 - σ * κ[i]) + B1 * ω == 0)
    @NLconstraint(m, ϕ2PDE[i=1:n], ϕ2pp[i] * σX[i]^2 / 2 + ϕ2p[i] * (μX[i] + σX[i] * (σ - κ[i])) + ϕ2[i] * (μ - r[i] - δ - δ2 - σ * κ[i]) + B2 * ω == 0)

    # Objective
    @NLobjective(m, Min, 1.)
    solve(m)