Like this
%Mod_Euler_method Modified Euler's method
% [t, w, h] = euler_method(f, a, b, alpha, n) performs Modified Euler's method for
% solving the IVP y' = f(t,y) with initial condition y(a) = alpha
% taking n steps from t = a to t = b.
j_ht = 1.75; % Jumper height in Metres
H = 74; % Height of jump point (Metres)
m = 80; % Mass of jumper (Kilograms)
g = 9.8; % Gravitational acceleration (m/s^2)
c = 0.9; % Drag coefficient (Kilogram/Metres)
k = 90; % Spring constant of bungee rope (Newton/Metres)
n = 10000; % Number of iterations
a = 0; % Start time
b = 60; % End time
D = 31; %Height of bridge (Metres)
C = c/m;
K = k/m;
L = 25; %length of rope in metres
h = (b-a)/n; % calculate h
t = a:h:b; % create time array t
% Create function handles for the system of ODE's
fv = @(v, y) v; % You don't need t here as you don't use t in the Euler iterations
fy = @(v, y) g - C*abs(v)*v - max(0, K*(y - L)); % ditto
% Also, make sure your arguments are in the same order here as when you
% call these functios in the modified Euler routine below!
y = zeros(1, n+1); % initialise y array
v = zeros(1, n+1); % initialise v array
% perform modified euler iterations
for i = 1:n
yk1 = h*fv(v(i), y(i));
vk1 = h*fy(v(i),y(i));
yk2 = h*fv(v(i) +yk1, y(i) + vk1);
vk2 = h*fy(v(i) +yk1, y(i) + vk1);
y(i+1) = y(i) + 1/2*(yk1+yk2);
v(i+1) = v(i) + 1/2 * (vk1 + vk2);
end
plot(t, y),grid
