Puzzling result from boundary condition code in Julia BVP solver

I am trying to solve a boundary value problem in Julia, following the example found here, using the BoundaryValueDiffEq package. In the boundary condition function, the example requires a for loop to update each index individually, à la

function bc1!(residual, u, p, t)
    for i in 1:n
        residual[i] = u[end][i] - 10
    end
end

I would like to use the following code, which should be more efficient:

function bc1!(residual, u, p, t)
    residual = u[end] .- 10
end

Though the resulting value of residual is the same for both versions of the code, the solver gives the correct result in the first case and an incorrect result in the second case.

All I can think of is that there is some difference between updating residual index by index and assigning a new vector to it, even if the result is identical in value and in type. Why is this the case, and is it possible to make the code more efficient while preserving the correct result?

Here is the full code in case it helps.

using BoundaryValueDiffEq, Plots

n = 3
f(t) = .1
F(t) = .1*t

function du!(du,u,p,t)
    fn(i) = 1/(u[i]-t)
    for i in 1:n
        du[i] = 1/(n-1)*F(u[i])/f(u[i])*((2-n)/(u[i]-t)+sum(map(fn,
        vcat(1:i-1,i+1:n))))
    end
end

function bc1!(residual, u, p, t)
    #residual = u[end] .- 10
    for i in 1:n
        residual[i] = u[end][i]-10
    end
end

# exact solution
xvals = LinRange(0,20/3,200)
yvals = 1.5*xvals

# solving BVP
tspan = (0.0,20/3)
bvp1 = BVProblem(du!, bc1!, 10*ones(Int8,n), tspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt=.2)

# plotting computed solution vs actual solution
plot(sol1,vars=(0,1))
plot!(xvals,yvals,label="Exact solution")

You overrode the array instead of mutating it. You need to use .= to update it in-place.

function bc1!(residual, u, p, t)
    residual .= u[end] .- 10
end

or safer:

function bc1!(residual, u, p, t)
    @. residual = u[end] .- 10
end