%Euler Method using 3 fns
%1) Driver, 2) RHS, 3) Euler's method
%  Andy Ruina,  Feb 15, 2021
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set up, driver
% parameters, ICs, time span, call the solver
t_end = 2*pi;  npoints = 20001;
tarray = linspace(0,t_end,npoints);
x0 = 1; v0 = 0;  z0 = [x0;v0];
p.k = 1;  p.m = 1; 

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Euler's method
% Appy the rule znew = znow + zdot * h
function [tarray zarray] = myEulerSolver(eqn, tarray, z0, p);
znow = z0;
npoints = length(tarray);

zarray = zeros(npoints, length(z0));

zarray(1,:) = znow'; 

for i = 1:npoints-1
  tnow = tarray(i);
  zdot = eqn(tnow, znow,p);
  h = tarray(i+1) - tarray(i);
    
  znew = znow + zdot*h; %The key thing in this whole process 

  zarray(i+1,:) = znew';
  znow = znew;

end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RHS function
% Our differential eqn in first order form
function  zdot =  myrhs(t, z, p)
    x= z(1); v= z(2);
    k = p.k; m= p.m;
    Force  = -k*x ;
xdot =     v;
vdot = Force/m;

zdot = [xdot;vdot];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%