This is the caller code. You can call it 'KS_monotonic_caller.m'. Copy the below script and paste the code in and save it in an appropriate directory.

% KS_MONOTONIC_CALLER will solve the K-S equation for the case of the
% monotonic shock conditions
% -------------------------------------------------------------------
%
%   Written 19 Mar 2021 for the INI Spring School
%           EXPONENTIAL ASYMPTOTICS FOR PHYSICAL APPLICATIONS

clear

ep = 0.05;                   % Set epsilon value
zmin = -25; zmax = 12;      % Set domain

% Define the initial condition at infinity 
A = 2;
ubc = @(z) 1 - A*exp(-2*z);
upbc = @(z) -2*(-A)*exp(-2*z);
uppbc = @(z) 4*(-A)*exp(-2*z);
ic = [ubc(zmax); upbc(zmax); uppbc(zmax)];

% Define the differential equation
fwd = @(t,Y) KSode(t,Y,ep);

% Solve the ODE from z = zmax going backwards to z = zmin
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-9);
[z, Y] = ode45(fwd, [zmax, zmin], ic, options);
u = Y(:,1);

figure(1)
hold all
plot(z, u);
plot(z, tanh(z), 'k--');
xlabel('z'); ylabel('u(z)');
ylim([-5,5]);
title('Monotonic shock solution (u0 = tanh(z) shown dashed)');
drawnow

The caller code will require 'KSode.m'. Save it and run it from the same directory. Make sure you name it with the correct upper and lower-case characters.

function Yp = KSode(z, Y, ep)
% KSODE provides the first-order differential equation definition for
%           ep^2 u''' + (1 - 4 ep^2) u' = 1 - u^2

u = Y(1);
up = Y(2);
upp = Y(3);

uppp = (1 - u^2 - (1 - 4*ep^2)*up)/ep^2;

Yp = [up; upp; uppp];

Again, you can copy and paste the below script into a Matlab file and call it, e.g. 'KS_osc_caller.m'. It will make use of the 'KSode.m' code above.

% KS_OSC_CALLER will solve the K-S equation for the case of the
% oscillatory shock conditions
% -------------------------------------------------------------------
%
%   Written 19 Mar 2021 for the INI Spring School
%           EXPONENTIAL ASYMPTOTICS FOR PHYSICAL APPLICATIONS

z0 = 1.24;
phi = 2.60586555;
% phi = 2.59;
% phi = 2.5;

ep = 0.05;
zmin = -6;
zmax = 12;

gam = sqrt(1/ep^2 - 1);
bc = @(z) -1 + exp(-(z-z0))*sin(gam*(z-z0)+phi);

bcp = @(z) -exp(-(z-z0))*sin(gam*(z-z0)+phi) + ...
    gam*exp(-(z-z0))*cos(gam*(z-z0)+phi);

bcpp = @(z) exp(-(z-z0))*sin(gam*(z-z0)+phi) ...
    - 2*gam*exp(-(z-z0))*cos(gam*(z-z0)+phi) ...
    -gam^2*exp(-(z-z0))*sin(gam*(z-z0)+phi);
ic = [bc(zmax); ...
    bcp(zmax); ...
    bcpp(zmax)];


fwd = @(t,Y) KSode(t,Y,ep);

mytol = 1e-10;
options = odeset('RelTol', mytol, 'AbsTol', mytol);
sol = ode113(fwd, [zmax, zmin], ic, options);

figure(1)
hold all
plot(sol.x, sol.y(1,:));
ylim([-5,5]);
drawnow