I am currently trying to solve the following problem in matlab- I am modelling a system of equations which I beleive follow this system of equations: $$x''=-GM\frac{x}{x^2+y^2}-\frac{\rho C_d A}{2m}x'\sqrt{x'^2+y'^2}$$ $$y''=-GM\frac{y}{x^2+y^2}-\frac{\rho C_d A}{2m}y'\sqrt{x'^2+y'^2}$$

This has the initial conditions and constants:

$x(0)=10000, y(0)=\alpha+R,x'(0)=7000,y'(0)=0,m=2000,A=12,C_d=1.5,R=2574.73\times10^3,M=1.3452\times10^{23},G=6.67\times10^{-11}$

This is made harder by the fact that $\rho$ is not a constant and in fact varies based on $(x,y)$

The aim of the project is to use the Runge-Kutta 4 method to solve these equations, so my first thought is to change this into a system with four coupled equations of single order ODEs. That went as follows:

$$z_1=x,\,z_1'=z_2$$ $$z_2=x',\,z_2'=-GM\frac{z_1}{z_1^2+z_3^2}-\frac{\rho C_d A}{2m}z_2\sqrt{z_2^2+z_4^2}$$ $$z_3=y,\,z_3'=z_4$$ $$z_4=y',\,z_4'=-GM\frac{z_3}{z_1^2+z_3^2}-\frac{\rho C_d A}{2m}z_4\sqrt{z_2^2+z_4^2}$$ $$z_1(0)=10000,z_2(0)=7000,z_3(0)=\alpha+R,z_4(0)=0$$

Any suggestions would be great, thanks in advance!


I have so far defined a function in matlab:

function dz = stateDeriv(t,z) 
  dz1 = z2*dt 


  dz = [dz1; dz2; dz3; dz4];

And we are given that:

function [rho] = profileTitan(z) 
%% profileTitan  Atmospheric density profile model for Titan  
% [RHO] = profileTitan(Z) outputs a vector RHO of modelled density values  
% at the corresponding altitude values in the input vector Z. Input  
% altitudes must be specified in meters and output densities are given in 
% kg/m^3. 
% This function uses the model of: 
% R.V. Yelle et al., "Engineering models for Titan's atmosphere", In Proc.  
% Hugyens Science, Payload, and Mission, European Space Agency Conference,  
% ESA SP-1177, pp. 243-256, Noordwijk, The Netherlands, 1997. 
%% Atmospheric model data for Titan
% Altitude samples 
z0 = [...
     1e-6 ...
     0.5 ...
     1 ...
     2 ...
     5 ...
     10 ...
     20 ...
     50 ...
     100 ...    
     200 ...
     500 ...
     1000 ...
     1300 ...
     ] * 1e3; % m

% Modelled density values from (R.V. Yelle et al., 1997) 
rho0 = [...
     0.579e-2 ...
     0.569e-2 ...
     0.559e-2 ...
     0.539e-2 ...
     0.478e-2 ...
     0.385e-2 ...
     0.237e-2 ...
     0.414e-3 ...
     0.296e-4 ...
     0.234e-5 ...
     0.101e-7 ...
     0.409e-11 ...
     0.168e-12 ...
     ]; % g/cm^3

% Convert density from g/cm^3 to kg/m^3 
rho0 = rho0 / 1e3 / ((1e-2)^3); 

%% Interpolate data at requested altitude(s)
% Interpolate data logarithmically for accuracy 
z0log = log(z0);
rho0log = log(rho0);

zlog = log(z);
% Handle negative altitudes zlog(imag(zlog)>0) = 0;
% Interpolate 
rholog = interp1(z0log,rho0log,zlog,'linear',-Inf); 
rho = exp(rholog);

% Assign infinite density for negative altitudes 
I = find(z<0); rho(I) = Inf;

For context, I have a spaceship that is aiming to start an eliptical orbit around a planet, centred at $(0,0)$. The forces we have to take into account are due to gravity and drag, since there is no thrust during this model. For this I used the formulas: $$F_d=\frac{1}{2}\rho v^2 C_d A$$ $$F_g=\frac{GMm}{r^2}$$ now to convert these into a vector form, so I can then split it up into two different equations, I defined: $$p=\begin{pmatrix} x\\y \end{pmatrix}$$ since drag will act in the opposite direction to velocity and gravity will act in the opposite direction to position, we will get: $$\vec{F_d}=-\frac{1}{2}\rho C_d A |\dot{p}|^2 \hat{\dot{p}}=-\frac{1}{2}\rho C_d A |\dot{p}|\vec{p}$$ $$\vec{F_g}=-\frac{GMm}{|p|^2}\hat{p}=-\frac{GMm}{|p|^3}\vec{p}$$ $$\Sigma F=m\ddot{p}$$

  • $\begingroup$ What exactly is your problem now? You did all that was necessary, now code it into one function that computes $z'$ from $z$ and insert it into a vector-enabled RK4 implementation of your choice. $\endgroup$ Commented Nov 11, 2019 at 15:05
  • $\begingroup$ My problem is that I have only seen the Runge Kutta used to solve with a single variable initially (transformed into two first orders, z1 and z2). $\endgroup$
    – Henry Lee
    Commented Nov 11, 2019 at 15:18
  • $\begingroup$ Runge Kutta doesn't do anything different with systems vs single equations, so the only thing wrong with what you wrote there is that you wrote, for example, z1 instead of z(1). Matlab will complain that z1 isn't defined. $\endgroup$
    – Ian
    Commented Nov 11, 2019 at 15:25
  • $\begingroup$ @Ian z1,z2,z3,z4 are newly defined variables, rather than different points $\endgroup$
    – Henry Lee
    Commented Nov 11, 2019 at 15:27
  • $\begingroup$ @HenryLee In your function declaration those aren't going to be defined in the scope you use them, unless they are global. dz1 etc. are fine, since those are being defined there. $\endgroup$
    – Ian
    Commented Nov 11, 2019 at 15:35

1 Answer 1


As @Ian commented, you need to follow the matlab array semantics. Additionally, the step size is not part of the derivative function. Using intermediary variables for the building blocks of the equation you can implement this as

function dz = stateDeriv(t,z) 
  r = sqrt(z(1)^2+z(3)^2);
  [rho] = profileTitan([r-R]);
  gravity = G*M/r^3;
  friction = (rho*Cd*A)/(2*m)*sqrt(z(2)^2+z(4)^2)
  dz(1) = z(2)
  dz(2) = -gravity*z(1)-friction*z(2)

  dz(3) = z(4)
  dz(4) = -gravity*z(3)-friction*z(4)

Then use this with the standard ode45 or ode4 from this video tutorial (see also How do you utilise the Runge Kutta method to calculate displacement and velocity from acceleration?)


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .