since you're dealing with a BVP, it's not a good choice to reduce to a first order system. That's because you should use a finite difference approach. Given a uniform grid from $0$ to $1$, with $N$ equally-spaced knots, i.e. $x_i=x_{i-1}+nh$, where $h$ is the discretization step, you can discretize the second and the first derivative (assuming knowledge of enough regularity of your solution) with finite differences, i.e.:
$u'(x_i) \approx \frac{u_{i+1}-u_{i-1}}{2h}$
$u''(x_i)\approx \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$
for $2 \leq i \leq N-1$, with the "convenction" that $u_i \approx u(x_i)$.
If you don't care for a moment about the boundary conditions, you can easily see (just put some values for $i$) that you can build a tridiagonal matrix which discretize the second derivative: so with a matrix-vector multiplication, you can get an approximation of $u''(x_i)$ for every $i$. Similarly for the first derivative. The entries will be made by 1,-2,1, for $u''(x_i)$, and -1,1 for $u'(x_i)$.
You can look here to see how they are built.
In a MatLab enviroment, you can build them with the commands
A = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),N,1));
B = toeplitz(sparse(2,1,-1/(2*h),m,1), sparse(2,1,1/(2*h),N,1));
(I used sparse because you should use a sufficient number of knots), but it's not important in order to understand how the matrices are)
In this way the original problem, after the above discretization, can be rewritten as
$A \vec{u} + \vec{u} (B \vec{u}-\vec{1})=0$
where $A,B$ are the matrices defined above (A for second derivative, B for first derivative). Of course the above equation is not much formal since I should give a meaning to what $\vec{u} (B-\vec{1})$ is. But if you look component-wise the discretization, you can see that it's just a component-wise multiplication. So, the correct way to write the discretized problem is
$A \vec{u} + diag(\vec{u}) (B \vec{u} -\vec{1})=0$
where you see immediately that the dimensions are OK.
Assuming you have imposed in the right way the boundary conditions (it's quite easy: just change first and last rows), in order to find $\vec{u}$ you have to solve a system of non-linear equations defined as $F(u)=A \vec{u} + diag(\vec{u}) (B \vec{u}-\vec{1})=0$.
The solutions $\bar{u}$ of $F(u)=0$ will be the numerical solution of you BVP
As you stated in your answer, Newton's method is the right choice, but we have to compute the Jacobian matrix of $F: \mathbb{R}^N \rightarrow \mathbb{R}^N$, defined as
$(JF)_{ij}=\frac{\partial F_i}{\partial u_j}$
The Jacobian of $A\vec{u}$ is just $A$ (think about the scalar case).
The term $\vec u$ has as Jacobian the identity matrix $I_N$ of dimension $N$ (differentiate component-wise the vector $(u_1,\ldots,u_n)$ w.r.t $u_1$ to $u_n$)
The last (and more difficult) term is $u u'$, where we have $diag(\vec u) (B\vec{u})$
Here it's more easy to look at what is the i-th component of this vector:
$u_i \cdot (\frac{u_{i+1}-u_{i-1}}{2h})$
In order to compute the Jacobian, we need to differentiate each component wrt the variables $u_1$ to $u_n$. But it's easy to see that here the only non vanishing terms are the ones which corresponds to positions $i-1,i,i+1$, and then the resulting matrix is tridiagonal.
More precisely we have this entries $\frac{-u_{i}}{2h}, \frac{u_{i+1}-u_{i-1}}{2h}, \frac{u_{i}}{2h}$. I just differentiate the above expression w.r.t $u_{i-1},u_i,u_{i+1}$
Since I like to build matrices, the resulting tridiagonal matrix can be written as $diag(B \cdot \vec{u})+diag(\vec{u})\cdot B$ (you can check).
Putting together all the terms, we have the following Jacobian matrix
$JF(\vec{u})=A-I_N+diag(B \cdot \vec{u})+diag(\vec{u})\cdot B$
and we can perform our Newton's method with a suitable initial guess $u_0$.
Of course, boundary conditions are really important and need to be imposed before the application of Newton's method.
EDIT [28/3/19]
The following runable Octave code solves your problem. Actually, I set $a=4,b=1,\delta=0.05$, but you only need to replace them with their effective values.
m=200; %grid points
h=1/(m-1);
x=linspace(0,1,m)';
%parameters corresponding to a,b,delta
a=4;
b=1;
delta=0.05;
A = toeplitz(sparse([1,2],[1,1],[-2,1]/(h^2),m,1));
B = toeplitz(sparse(2,1,-1/(2*h),m,1), sparse(2,1,1/(2*h),m,1));
F=@(u) [u(1)-a;(delta*A*u+diag(u)*(B*u-1))(2:m-1);u(m)-b];
JF=@(u) [[1,zeros(1,m-1)];(delta*A+diag(B*u)+diag(u)*B)(2:m-1,1:m);[zeros(1,m-1),1]];
%Newton's method
u0=linspace(a,b,m)'; %initial guess
tol=h^2/100; %second order FD
res= -JF(u0)\F(u0); %residue
while(norm(res,inf)>tol)
u0+=res;
res=-JF(u0)\F(u0); %main loop of Newton's method
endwhile
u0+=res;
plot(x,u0,'*')
I also checked the right order of convergence, you can see it in the following plot

I hope I have been helpful.