Consider a Boundary Value Problem:
$$ \delta u''+u(u'-1) =0 \Leftrightarrow u''=\frac{-u(u'-1)}{\delta}=:f(t,u',u), \\ u(0)=a, u(1)=b $$
$\delta,a,b$ are known parameters. I want to implement Newton's method on a uniform grid. How do I derive the Jacobian?
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.
$$G(u)=\delta u''+u(u'-1) =0 \\ u(0)=a, u(1)=b$$ If you wish to calculate implicitly the Jacobian before applying any discretization scheme on your PDE first (e.g. so you can reuse the code of a linear solver already available to you), you need the linearization of the nonlinear operator $G(u)$. This is the Frechet derivative of $G$, defined as $$G_u \Delta\equiv\lim_{\epsilon\rightarrow0}\frac{G(u+\epsilon \Delta)-G(u)}{\epsilon}$$ Here, $\Delta$ is also a function.
Its derivation is straight forward, plug in $G(u)$ into the definition of the Frechet derivative and take the limit in $\epsilon$. This gives you $G_u \Delta = \delta \Delta''+ u \Delta' + (u'-1) \Delta$, now a linear operator in $\Delta$.
In the Newton-Kantorovich iteration one uses the generalized Taylor series $G(u+\Delta)=G(u)+G_u\Delta+O(\Delta^2)$ at a solution $u^{(i)}$ of some previous iteration step which solves already the boundary conditions, and solves a linear PDE for $\Delta$ with homogeneous boundary conditions $\Delta(0)=\Delta(1)=0$, $$0=G^{(i)}+G_u^{(i)} \Delta$$ Once you have found the solution for $\Delta$ you update $u^{(i)}$ to $$u^{(i+1)}=u^{(i)}+\Delta$$.
This is repeated until $\Delta$ is small in some sense, e.g. $\int_\Omega\Delta^2 <\varepsilon$.
So the Newton part is more or less identical to VoB answer. We simply skip the all the error prone component wise calculations and go right to the correct definition of the derivative, so we can reuse any well tested well written linear PDE solver.
