How do I define the boundary condition function for bvp4c for shock wave on a cone?

I am trying to solve the shock wave problem on a cone (see the methodology and figure with variables attached). The problem is driven by an ODE with boundary conditions, so I am relying on the Matlab function bpv4c. I am working with this tutorial and I try to adapt it to my problem.

I am good with the definition of odefun (equations (1), (2) and (3)), but I don't see at all how to code the boundary condition function bcfun (equations (4), applicable at Theta_s and (5), applicable at Theta) and the initial solution solinit.

Your help would be greatly appreciated.

Input parameters

- Theta_s: typically between 5 ° and 20 °

- Theta: current angle, greater than Theta_s

- Gamma: typically 1.4

- u_s: between 0 and 1

Output parameters

- u and v, for angles between Theta_s and Theta

ODE Function

function dV_dTheta = shock_taylor_maccoll(Theta,V,Gamma)

dV_dTheta = NaN(2,1); % Pre-allocation of derivative vector

A = (Gamma - 1) / 2 * (1 - V(1) ^ 2 - V(2) ^ 2); % Intermediate calculation

dV_dTheta(1) = V(2); % Derivative of radial velocity equals angular velocity
dV_dTheta(2) = (A * (V(1) + V(2) * cotd(Theta))) / ...
    (V(2) ^ 2 - A) - V(1); % Derivative of angular velocity
end

[Christian Soulard]

Update with the right figure for the shock on the cone, sorry for the mistake.