MateLab + Latex

function [t, y] = solve_system_euler(f, y_start, t_start, t_final, dt) % solves a general system of first-order ODEs % y1' = f1(t, y1, y2, ... , ym) % y2' = f2(t, y1, y2, ... , ym) % ... % ym' = fm(t, y1, y2, ... , ym) % with initial conditions % y1(t_start) = y_start_1, % y2(t_start) = y_start_2, % ... % ym(t_start) = y_start_m, % using the Euler Method % % input: % f = right-hand side (vector-valued function) % y_start = initial value (vector of initial conditions) % t_start = starting time % t_final = final time % dt = desired time step % output: % t = array of times % y = matrix containing numerical solution (each column corresponds to % solution at certain time t_i) % i.e., the output has the following form % t = [ t_1, t_2, t_3, ... , t_n ] % y = [ y1_1, y1_2, y1_3, ... , y1_n; % y2_1, y2_2, y2_3, ... , y2_n; % ... , ... , ... , ... , ... ; % ym_1, ym_2, ym_3, ... , ym_n ] % % where n is the total number of steps in time % and m is the size of a system % Note that f has to be a vector-valued function of t and y, % where y is a vector % compute total of number of iterations num_steps = ceil((t_final-t_start)/dt); % determine size of the system m = length(y_start); % allocate arrays y = zeros(m, num_steps); t = zeros(1, num_steps); % impose initial conditions t(1) = t_start; y(:,1) = y_start; % symbol ':' is used to refer to the entire column % this is equivalent to: % for i = 1:m % y(i,1) = y_start(i) % end n = 1; while t(n) < t_final if t(n) + dt > t_final dt = t_final - t(n); end t(n+1) = t(n) + dt; y(:, n+1) = y(:, n) + dt * f(t(n), y(:, n)); % this is equivalent to: % % slopes = f(t(n), y(:,n)); % for i = 1:m % y(i,n+1) = y(i,n) + dt * slopes(i); % end n = n+1; end end