Differential equation concerning echolocation

I am in the process of trying to solve a differential equation using MATLAB.

Given information is the following: $$c(z) = 4800 - 20.2090 + \frac{17.3368z}{1000}+272.9057e^\frac{-0.7528z}{1000}$$ $$\frac{d^2z}{dx^2} = -q_0\frac{c'(z)}{c(z)^3}$$ $$q_0=(\frac{c(z_0)}{cos(B_0)})^2$$ $$x = 0, z = z_0, \frac{dz}{dx} = tan(B_0)$$ B0 is the angle between the line z = 2000 and the ray in the start point.

I have solved the equation using MATLAB's ode45 command, but now I'm supposed to use the secant method to find which initial values of $B_0$ also give me an end value of $z = 2500$ at $x = 151900$. I've been told that there's 4, since I'm supposed to use the $\Delta t$ between the arrivals to pinpoint the initial position, but no matter how I try I cannot figure out how to use the secant method in this particular instance.

Currently my program looks something like this:

angle = (-10:0.0035:14);

z_0 = 2000;

error_tolerance = 1;    % Väljer hur långt ljudvågen
it = 1;

figure
for i = 1:length(angle)

B_0 = tand(angle(i));
q_0 = (c_sound(z_0)/cosd(angle(i)))^2;
[X,Z] = ode45(@soundwave_function, [0 distance], [z_0 B_0]);
if abs(Z(end,1)-2500) <= error_tolerance

saved_angles(it,1) = angle(i);

it = it+1;

end

end

angles_str = num2str(saved_angles);

global saved_angles

if numel(saved_angles)>1

for i = 1:numel(saved_angles)
disp([num2str(saved_angles(i)), ' degrees'])
end

else
end


This does find the angles where $z(151900) \approx 2500$ in the interval specified, -10 to 14 degrees (with 0 degrees counting as z = 2000).

I don't know how to approach the problem using the secant method, and since I'm supposed to use that solution for the second bit I don't really know how to continue.

Is there anyone out there that knows how to apply the secant method in this specific case? Would be much appreciated!

You need to use a shooting method. Let's abstract it a bit.

Let's define $F(x_{ic})$ to be a function that takes your initial conditions as its argument and outputs the value of your solution at $x = 151900$.

In other words, $F$ is a solution trajectory of your ODE, and we might think of $F$ as a single run of ode45, executed with initial conditions $x_{ic}$.

If we want to compute $F(x_{ic}) - 2500 = 0$, we can find initial values that satisfy the problem. Therefore, we simply want to write a root finding algorithm.

Newton's method is a good fit, but it is hard to compute $\frac{dF}{dx_{ic}}$. Instead, we might use the secant method, also known as Broyden's method.

It is not so hard: simply implement a function for Broyden's method, and wrap your ODE solver in a function that takes initial conditions as its arguments. For instance:

function z151900 = myODE(IC)
[x,z] = ode45(myode,[0 151900],IC);
z151900 = z(end);
end;

function xic = SecantMethod(F,IC1,IC2)
% your secant method iteration using two initial guesses
end;

IC1 = ic1;
IC2 = ic2;
SecantMethod(z151900,ic1,ic2);


I can post a more concrete example later, or you can search my answer history for "Broyden's method".