Matlab code using ode45 difficulty

[zfolwick]

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.

 

[Nylex]

I've gotten nothing but errors.

Don't you think it would be helpful if you post the errors?

 

[GreenPrint]

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?

 

[Nylex]

It also wouldn't hurt to use

 tags to preserve the indentation, to help readability.