%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function prel1q3()        
 tspan=[0 6.2];   n = 100;   
 z0   = [1 0]';   p =0;
 [tarray zarray] = eulmeth(tspan,z0,n,p);
 u = zarray(:,1);  
 w = zarray(:,2);
 plot(u,w,'r', tarray, u,'b');  
 axis('equal')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tmat zmat] = eulmeth(tspan,z0,n,p)
tmat   = linspace(tspan(1),tspan(2),n+1); 
h      = tmat(2)-tmat(1);
zmat   = zeros(n+1,length(z0)); 
zmat(1,:) = z0';
 for i=1:n;
    z           = zmat(i,:)'; 
    t           = tmat(i);
    znew        = z + h*rhs(t,z,p); 
    zmat(i+1,:) = znew';
 end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  zdot = rhs(t,z,p)
  z1   = z(1);  
  z2   = z(2);
  zdot = [ z2;   -z1 - p*z2 ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%