```
function BVPattempt
r=10; % radius of circular flow (kpc)
alpha=0.11667; % angle of spiral inclination (radians)
solinit = bvpinit(linspace(0,alpha*r*pi,1000),[13.4171 115]); % not sure about the last argument, its a constant guess for the solution across the grid
sol = bvp5c(@f,@ex1bc,solinit);
% The solution at the mesh points
x = sol.x;
y = sol.y;
figure;
plot(x,y(1,:)')
hold on
plot(x,y(2,:)')
legend('u','v');
xlabel('eta');
ylabel('u/v');
title('None');
%---------------------------------------------------------------
function dudy = f(y,u)
dudy = zeros(2,1);
% define physical constants
omega=25; % angular velocity at omega(r) (km /sec /kpc)
kappa=31.3; % epicyclic frequency (km /sec /kpc)
omegapattern=13.5; % angular pattern velocity (km /sec /kpc)
c=8.56; % speed of sound (km /sec)
r=10; % radius of circular flow (kpc)
alpha=0.11667; % angle of spiral inclination (radians)
A=72.92; % amplitude (km /sec)^2
v0=r*(omega - omegapattern);
u0=alpha*r*(omega - omegapattern);
dudy(1) = (u(1)/((u(1)^2) - (c^2)))*(2*omega*(u(2) - v0) + ((2*A)/(alpha*r))*sin(((2*y)/(alpha*r))));
dudy(2) = -((kappa^2)/(2*omega))*((u(1) - u0)/u(1));
%---------------------------------------------------------------
% Not sure how to impose two boundary conditions
function res = ex1bc(ya,yb)
%res = [ ya(1) - 13.4171
%yb(1) - 13.4171];
res = [ ya(1) - 13.4171
yb(1) - 13.4171];
```