%Euler Method using 3 fns
%Minimal BAD solution.  Andy Ruina,  Feb 15, 2021
% Minimized, not for clarity, but to show how little is really here.
% Do not imitate this!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tarray = linspace(0,2*pi, 101);
z0 = [1;0];
p.k = 1;  p.m = 1; 

[tarray zarray] = myEulerSolver(@myrhs, tarray, z0, p);
xarray  =   zarray(:,1); varray = zarray(:,2);
plot(xarray,varray); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Euler's method
function [tarray zarray] = myEulerSolver(eqn, tarray, z0, p);
z           = z0;  %z is present value of z
npoints     = length(tarray);
zarray      = zeros(npoints, length(z0));
zarray(1,:) = z'; 

for i = 1:npoints-1
  t             = tarray(i);
  zdot          = eqn(t, z,p);
  h             = tarray(i+1) - tarray(i);
  znew          = z + zdot*h; %The key thing in this whole process 
  zarray(i+1,:) = z';
  z             = znew;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RHS function
% Our differential eqn in first order form
function  zdot =  myrhs(t, z, p)
x= z(1); v= z(2);
xdot =     v;
vdot = -p.k*x/p.m;
zdot = [xdot;vdot];
%The 4 lines above could be replaced with a single line, but it is hard
%to read and bad practice.  Don't do this:
% zdot = [z(2);  -p.k*z(1)/p.m];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%