# Matlab code using ode45 difficulty

## Homework Statement

I've been trying to get this code to work:
ode45(@(t,y) fallode(t,y,B), [0, tmax], [31330, 0])

## Homework Equations

function dy = fallode(t, y, B)
% ODE function to model a fall over a large range of
% altitudes, reaching up to high subsonic Mach numbers.
% y(1) is altitude, y(2) is vertical velocity (positive up).
% Variation of gravity with altitude is ignored.
g = 9.8; % Earth gravity, m/s^2
R = 287; % Specific gas constant of air, J/kg*K
gamma = 1.4; % Ratio of specific heats of air, dimensionless
[T, rho] = atmos(y(1));
dy(1,1) = y(2);
dy(2,1) = -g + rho*B*y(2)^2/(2*sqrt(1 - y(2)^2/(gamma*R*T)));

--------------

function [T, rho] = atmos(ha)
% Calculate the temperature and density at a given absolute
% altitude ha in the US standard atmosphere model (lowest
% three layers only)
re = 6378.e3; % Radius of Earth, mean equatorial, meters
% Geopotential altitude
h = re*ha/(re + ha);
R = 287; % Specific gas constant for air, J/(kg*K)
g0 = 9.8; % Earth surface gravity, m/s^2
hbreak1 = 11e3; % Altitude of break between 1st and 2nd layers in model
hbreak2 = 25e3; % Altitude of break between 2nd and 3rd layers in model
for ii = 1:length(ha)
if (h <= hbreak1)
a = -6.5e-3; % Temperature lapse rate in troposphere, K/m
rho1 = 1.225; % Surface air density in standard atm. model, kg/m^3
T1 = 288.16; % Surface temperature in standard atm. model, K
T(ii) = T1 + a*h(ii); % Temperature at altitude, K
rho(ii) = rho1*(T(ii)/T1)^(-g0/(a*R)-1); % Density at altitude, kg/m^3
elseif (h <= hbreak2)
T(ii) = 216.66; % Temperature between 11km and 25km altitude, K
rhos = 0.3642; % Density at 11km altitude, kg/m^3
rho(ii) = rhos*exp(-(h(ii) - hbreak1)*g0/(R*T)); % Density at altitude, kg/m^3
else
4a = 3e-3; % Temperature lapse rate in stratosphere, K/m
rho1 = 0.0401; % Density at 25km altitude, kg/m^3
T1 = 216.66; % Temperature at 25 km altitude, K
T(ii) = T1 + a*(h(ii) - hbreak2); % Temperature at altitude, K
rho(ii) = rho1*(T(ii)/T1)^(-g0/(a*R)-1); % Density at altitude, kg/m^3
end
end

## The Attempt at a Solution

I've gotten nothing but errors. Any help troubleshooting would be appreciated.

Related Engineering and Comp Sci Homework Help News on Phys.org
I've gotten nothing but errors.
Don't you think it would be helpful if you post the errors?

Last edited:
Don't you think it would be helpful if you post the errors?
I concur, also why don't you just use fprintf to output to the display screen your results with units instead of using notes that don't get outputted to the display screen?

It also wouldn't hurt to use
Code:
 tags to preserve the indentation, to help readability.