% Demo 3b in BAD compact form
clear

% PART 1: General Setup
p.m    = 1;    p.k = 1;                   
tarray = linspace(0, 2*pi, 100);      
z0     = [1;0];      

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

% PART II: Num. method
function  [tarray, zarray]   = myEulerSolver(eqn, tarray, z0, p)
ntimes = length(tarray);     
zarray = zeros(ntimes, length(z0));  
zarray(1,:) = z0';    
 for i  = 2:ntimes    
  h    = tarray(i) - tarray(i-1);       
  zarray(i,:) = zarray(i-1,:)+ eqn(tarray(i), zarray(i-1,:)',p)' *h;
  end 
end

% PART III:  The ODEs
function zdot = myrhs(t,z,p)
 zdot= [z(2); -(p.k/p.m)* z(1)];
end 