Solving System of 1st Order ODEs with Euler's method

I have a problem that requres taking a second order ODE and decomposing it into 2 first-order ODEs, then approximating a solution using Euler's explicit method. I already did the decomposition:

Here are the resultant ODEs:

y1' = y2

y2' = [ (5000+80t-0.161y2^2)*(32.2/(3000-80t)) ]

As you can see I have one dependent varialble y, and one independent variable t.

My question is how to write an Euler function file with 2 equations. I have an Euler function file from a textbook that takes care of a single ODE, but I want to solve a system of coupled ODEs. Below is the Euler code from the textbook for a single ODE, but I don't even know where to start in terms of writing my own code.

% function file 
function [x, y] = odeEULER(ODE,a,b,h,yINI) 

x(l) = a; y(l) = yINI; 
N = (b - a)/h; 

for i = l:N 
    x(i + 1) = x(i) + h; 
    y(i + 1) = y(i) + ODE(x(i) ,y(i))*h; 
end 

end

Answer 1

[James Tursa]

The easiest way is to treat your y as 2-element column vectors instead of scalars. E.g.,

% function file 
function [x, y] = odeEULER(ODE,a,b,h,yINI) 

x(1) = a; 
y(:,1) = yINI; % yINI needs to be initial 2-element [y1;y1'] vector
N = (b - a)/h; 

for i = l:N 
    x(i + 1) = x(i) + h; 
    y(:,i + 1) = y(:,i) + ODE(x(i) ,y(:,i))*h; 
end 

end

You would of course have to modify your ODE( ) function to return a 2-element column vector instead of a scalar.