I know this doesn’t directly answer your question, but BSplineKit.jl allows Robin-type boundary conditions.

The first thing you would need to do is to convert your boundary conditions to homogeneous BCs, for instance via the change of variables v = u - 1. Your equation then becomes:

The interface for solving boundary value problems with BSplineKit is currently not exactly high level. I can explain later how this works, but right now one would do it as follows:

using BSplineKit
using LinearAlgebra
ts = range(0, 1; length = 101)
B = BSplineBasis(BSplineOrder(4), ts)
# Homogeneous boundary conditions
bc_left = Derivative(1)
bc_right = Derivative(0) + Derivative(1) # Robin-type BC
# Create recombined B-spline basis, such that each element of the basis individually satisfies
# the specified homogeneous BCs
R = RecombinedBSplineBasis(B, bc_left, bc_right)
M″ = galerkin_matrix(R, (Derivative(0), Derivative(2)))
M = galerkin_matrix(R, (Derivative(0), Derivative(0)))
A = M - M″
y = galerkin_projection(x -> -one(x), R) # project right-hand side
vcoefs = A \ y # solve linear system for basis coefficients
v = Spline(R, vcoefs) # finally, create spline with the solution (can be evaluated, plotted, differentiated, ...)
## Verification
v′ = Derivative() * v
v″ = Derivative() * v′
v′(0) # = 0
v(1) + v′(1) # ≈ 0
# Verify graphically that v″ = v + 1
using CairoMakie
fig = Figure()
ax = Axis(fig[1, 1]; xlabel = L"x")
lines!(ax, 0..1, x -> v(x) + 1; label = L"v(x) + 1")
lines!(ax, 0..1, x -> v″(x); linestyle = :dash, label = L"v″(x)")
axislegend(ax; position = (0, 1), labelsize = 20)
fig

@jipolanco: I will see if I can use that for my actual problem, although I would like to use ApproxFun in priority.

To Sevi: My problem is a linear one, thus I don’t want to use Newton, but I will give it a try.

@jishnub: Actually, I tried something like this, but It didn’t work (I am new to Julia so I don’t know the intricacies of these operators yet). I tried for instance

d = Domain(0…1)
L = -D^2-I
u = [Neumann(d)+(x>0.5)*Dirichlet(d);L] \ [[0,1],0]

but that didn’t give the right solution. How would you impose Robin boundary conditions using a sum of Dirichlet and Neumann operators?

The operator \ accepts a tolerance, so instead of A \ b, you may call it as (A, b, tolerance=1e-16) (or whatever tolerance you want).

If you want higher precision, perhaps you could try using a BigFloat or a BigInt domain? This should also automatically bump the tolerances used. Something like d = Domain(big(0)..big(1)).