A onelayer model of the largescale circulation in an ocean (the Gulf Stream) was proposed by Ierley and Ruehr in 1986 [1]: $$ u''' \lambda((u')^2uu'')u+1=0, ~~u(0) =0, ~~x\in [0,\infty) $$ for some real $\lambda$. Boundary conditions are given as either $u(0)=u'(0)=0$ (rigid or noslip) or $u(0)=u''(0)=0$ (slippery or stressfree). In both cases we require $u(\infty) = 1$.
In order to solve the problem in Chebfun we'll need to truncate the domain to something suitable, say $[0, X]$, i.e, to make use of the socalled domain truncation (see for instance our papers [3] and [4]).
An integral result for this problem is fairly useful. It reads $$ I = \int_{0}^{\infty }[(u'')^2  3\lambda uu'u'']dx=\frac{1}{2}, $$ and it is obtained multiplying the equation by $u'$, integrating by parts and enforcing the boundary conditions. This result is valid for both types of boundary conditions. Using this integral result we optimise the value of the length $X$ above, and find that the accuracy of the Chebfun result comes close to machine precision.
We can set up the chebop and solve the differential equation with only a few lines of code (see [2] for details).
tic X = 35; dom = [0, X]; lambda = 0.1; op = @(u) diff(u,3)  lambda*( diff(u,1)^2  u*diff(u,2) )  u + 1; lbc = @(u) [u; diff(u,2)]; % stressfree BC rbc = 1; N = chebop(op,dom,lbc,rbc); [u,info] = N\0;
Here is what the solution looks like.
plot([u diff(u) diff(u,2)]) axis([0 20 1 1.5]) xlabel('x'), legend('u','u''','u''''','location','southeast') title('Slippery or stressfree b. c.')
The residuals are small:
N_residual = norm(N(u)) % residual of diffl. eq. lbcu = lbc(u); lbc_residuals = [lbcu{1}(0) lbcu{2}(0)] % residuals of left BC rbc_residual = u(end)  rbc % residual of right BC
N_residual = 4.580588942258870e10 lbc_residuals = 1.0e09 * 0.000033013148830 0.133907788023016 rbc_residual = 9.769962616701378e15
The Newton iteration has converged quadratically:
semilogy(info.normDelta,'m*') ylim([1e16 1e+01]) xlabel('iteration') ylabel('norm of Newton update') %
The Chebyshev coefficients of the solution decrease rapidly:
plotcoeffs(u)
Finally, the integral $I$ comes out with a very small error:
I = sum(diff(u,2)^2  3*lambda*(u*diff(u)*diff(u,2))) I_error = abs(I1/2)
I = 0.499999999999915 I_error = 8.482103908136196e14
We solved the problem for several values of $X$ and found that the minimal error in $I$ occurs with $X\approx 35$.
total_time_for_this_example = toc
total_time_for_this_example = 3.629420000000000
References

G. R. Ierley and O. G. Ruehr, Analytic and numerical solutions of a nonlinear boundarylayer problem, Stud. Apl. Math. 75:136 (1986).

L. N. Trefethen, A. Birkisson, and T. A. Driscoll, Exploring ODEs, SIAM, 2018.

C. I. Gheorghiu, Pseudospectral solutions to some singular nonlinear BVPs, Numer. Algor. 68 (2015), 114, DOI: 10.1007/s110750149834z.

C. I. Gheorghiu, Spectral collocation solutions to systems of boundary layer type, Numer. Algor. 73 (2016), 114, DOI:10.1007/s1107501500836