MATLAB: Can anybody help me in figuring out this coding? It gives me the following error : (Undefined function ‘rfk45’ for input arguments of type ‘function_handle’. Error in Example_11_03 (line 37) [t,f] = rfk45(@rates,tspan,f0);. )

Question

"function Example_11_03
clear all; close all; clc
deg = pi/180;
g0= 9.81;
Re = 6378e3;
hscale = 7.5e3;
rho0 = 1.225;
diam = 196.85/12*0.3048;
A = pi/4*(diam)^2;
CD = 0.5;
m0 = 149912*.4536;
n = 15;
Isp = 390;
T = 933910;
mfinal = m0/n;
me = Thrust/Isp/g0;
mprop = m0-mfinal;
tburn = mprop/me;
hturn = 130;
t0 = 0;
tf = tburn;
tspan = [t0,tf];
v0 = 0;
gamma0 = 89.85*deg;
x0 = 0;
h0 =0;
vD0 = 0;
vG0 = 0;
f0 = [v0; gamma0; x0; h0; vD0; vG0];
[t,f] = rfk45(@rates,tspan,f0);
v = f(:,1)*1.e-3;
gamma = f(:,2)/deg;
x = f(:,3)*1.e-3;
h = f(:,4)*1.e-3;
vD = -f(:,5)*1.e-3;
vG = -f(:,6)*1.e-3;
for i = 1:length(t)
    Rho = rho0*exp(-h(i)*1000/hscale);
    q(i) = 1/2*Rho*v(i)^2;
end
output
return
function dydt = rates (t,y)
dfdt = zeros(6,1);
v = y(1);
gamma = y(2);
x = y(3);
h = y(4);
vD = y(5);
vG = y(6);
if t < tburn
    m = m0-me*t;
    T = Thrust;
else
    m = m0-me*tburn;
    T = 0;
end
g = g0/(1+h/Re)^2;
rho = rho0*exp(-h/hscale);
D = 1/2*rho*v^2*A*CD;
if h<= hturn
    gamma_dot = 0;
    v_dot = T/m - D/m - g;
    x_dot = 0;
    h_dot = v;
    vG_dot = -g;
else
    v_dot = T/m - D/m - g*sin(gamma);
    gamma_dot = -1/v*(g-v^2/(Re+h))*cos(gamma);
    x_dot = Re/(Re+h)*v*cos(gamma);
    h_dot = v*sin(gamma);
    vG_dot = -g*sin(gamma);
end
vD_dot = -D/m;
dydt(1) = v_dot;
dydt(2) = gamma_dot;
dydt(3) = x_dot;
dydt(4) = h_dot;
dydt(5) = vD_dot;
dydt(6) = vG_dot;
end
function output
fprintf('\n\n---------------------------------\n')
fprintf('\n Initial flight path angle = %10g deg',gamma0/deg)
fprintf('\n Pitchover altitude        = %10g m',hturn)
fprintf('\n Burn time                 = %10g s',tburn)
fprintf('\n Final speed               = %10g km/s',v(end))
fprintf('\n Final flight path angle   = %10g deg',gamma(end))
fprintf('\n Altitude                  = %10g km',h(end))
fprintf('\n Downrange distance        = %10g km',x(end))
fprintf('\n Drag loss                 = %10g km/s',vD(end));
fprintf('\n Gravity loss              = %10g km/s',vG(end))
fprintf('\n\n---------------------------------\n')
figure(1)
plot(x,h)
axis equal
xlabel('Downrange distance (km)')
ylabel('Altitude (km)')
axis([-inf, inf, 0, inf])
grid
figure(2)
subplot(2,1,1)
plot(h,v)
xlabel('Altitude (km)')
ylabel('Speed (km/s)')
axis([-inf, inf, -inf, inf])
grid
subplot(2,1,2)
plot(t,gamma)
xlabel('Time (s)')
ylabel('Flight path angle(deg)')
axis([-inf, inf, -inf, inf])
grid
figure(3)
plot(h,q)
xlabel('Altitude (km)')
ylabel('Dynamic pressure (N/m^2)')
axis([-inf, inf, -inf, inf])
grid
end
end"

Answer 1

rfk45 in not a standard MATLAB function. You will either need write it yourself or get it somewhere else... there are several versions on MATLAB FEX, such as this one:

https://www.mathworks.com/matlabcentral/fileexchange/2182-numerical-methods-using-matlab--3e?focused=5042708&tab=function