# ode45 usage in this case?

The background is that I'm solving a problem in numerical analysis and I seem to be halfway having solved the non-linear system and now getting to the part where I should apply solving differential equation.

https://math.stackexchange.com/questions/467164/how-to-proceed-solving-this-problem

The code so far is

z=[0:500:4000 5000:1000:12000];
data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
x0=[0 0 0 -1];
opt = optimset('MaxFunEvals',1000);
p=lsqnonlin(fun,x0,[],[],opt);
fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t);
tt=linspace(0,12000,1000);
plot(z,data,'r-',tt,fitf(tt),'b-');


And the system of equations is stated in the math question. Now I suppose that I should define a function something like

function f=fsys(x, z)
f= ...


But first I suppose that I must rewrite the system from second order to first order which is just math theory, right? So I first rewrite the system to first order and then I define the function in matlab and finally something like

xend=...; z0=...
[X,Z]=ode45(@fsys,[0, xend],z0); zs=Z(end), plot(X,Z,'--')


Can you help me define fsys, tell me how to proceed, confirm or reject that I'm on the right way? Any hints if you don't want to help me but hint at how to learn matlab well enough to solve this "pretty simple math"?

I think I can transform the system in the usual way:

u = z'
u' = z''


The system becomes

u'=-q0*c'(z)/((c(z)^3)


Since z is a function of x, u and u' are also functions of x and thus this system till become a first-ordet system but how do I program the solution in matlab? How do I use that u=tan(b0)?

## 1 Answer

I take $p4=1$ and i get the other coefficients as $p_{1}=5.66$, $p_{2}=14.58$, $p_{3}=257.4$. So $c(z)$ and $c'(z)$ will be as

$c(z)=4800+5.66+(14.58)\frac{z}{1000}+(257.4)exp(-\frac{z}{1000})$ and

$c'(z)=(1287exp(\frac{z}{1000}))/5000+729/50000$

Let $z'=u$ and $z''=u'$ than we may write following differential eq.

$$\frac{d^{2}z}{dx^{2}}=\left(\frac{c(2000)}{cos(7.8)}\right)^{2}\frac{c'(z)}{c(z)^{3}}, \qquad z(0)=2000,\qquad z'(0)=tan(7.8)$$

as

\begin{align*} z' & =u,\qquad z(0)=2000\\ u' & =\left(\frac{c(2000)}{cos(7.8)}\right)^{2}\frac{c'(z)}{c(z)^{3}},\qquad u(0)=tan(7.8) \end{align*}

So the codes for the problem should be as

    function dZ=sys(x,Z)
c=@(z)4800 + 5.66 + (14.58)*z/1000+ (257.4)*exp(-z/1000); % c(z)
c=c(2000);

% Z(1):=z
% Z(2):=u
dZ=zeros(2,1);    % a column vector
dZ(1)=Z(2);
dZ(2)=(c/cosd(7.8))^2*((1287*exp(Z(1)/1000))/5000 + 729/50000)/...
(4800 + 5.66 + (14.58)*Z(1)/1000+ (257.4)*exp(-Z(1)/1000))^3;
end


We can call the above function by the following codes and plot the solution

>> x=0:0.5:3000;
>> [X,Z]=ode45(@sys,x,[2000 tand(7.8)]);
>> plot(X,Z(:,1),'r',X,Z(:,2),'b')  %Z(:,1) is z(x) and Z(:,2) is z'(x).


You see just z below

• Thanks! But we've gotten different coefficients. I have to check that the coefficients are correct. The matlab I run is z=[0:500:4000 5000:1000:12000]; data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990]; fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data; x0=[0 0 0 -1]; opt = optimset('MaxFunEvals',1000); p=lsqnonlin(fun,x0,[],[],opt); fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t); tt=linspace(0,12000,1000); plot(z,data,'r-',tt,fitf(tt),'b-'); and then p(4) is not 1, is it something else? – Niklas R. Aug 31 '13 at 12:06
• The coefficients i got would be more precise actually i didn't focus on coeffs compeletely, because you were trouble with conversion of the original ode to the system of ode so i focused on the ode45 usage. Anyway you can use your coefficients also. Just replace your coeffs in the codes. – Ömer Aug 31 '13 at 21:33