function fig227
clc
clf
disp('Page 163 Fig2.27: harmonic oscillator System')
disp('dy/dt = v')
disp('dv/dt = -y')
disp('To be shown: vector fields (unnormalized)')
pause;

% Fig. 2.27
% Harmonic oscillator system
ymin = -1;
ymax = 1;
vmin = -1;
vmax = 1;
dy = 0.1;
dv = 0.1;
[y,v] = meshgrid(ymin:dy:ymax,vmin:dv:vmax);
u0 = v;
v0 = -y;
leng = sqrt(u0.^2+v0.^2);
%
clf
% unnormalized version
quiver(y,v,u0,v0,0.5)
axis equal tight
xlabel 'y'
ylabel 'v'
title 'Vector Field for Harmonic Oscillator System'
%
disp('To be shown: vector fields (normalized)')
pause
% normalize version
quiver(y,v,u0./leng,v0./leng,0.4)
axis equal tight
xlabel 'y'
ylabel 'v'
title 'Vector Field for Harmonic Oscillator System'
%
disp('To be shown: equilibrium points')
pause
% solve for the equilibrium points
[y0,v0] = solve('v0=0','-y0=0', 'y0','v0');
[y0 v0]
hold on;
plot(double(y0),double(v0),'go')
%
disp('To be shown: solution curves with initial condition (y,v)(0)=(0.3, 0)')
% solve for the solution curve
fode = inline('[y(2); -y(1)]', 't', 'y');
y10 = 0.3;
y20 = 0.0;
t0 = 0;
tmax = 20;
[t,ysol] = ode45(fode,[t0 tmax], [y10 y20]);
hold on
% solution curve in phase space
plot(ysol(:,1),ysol(:,2),'r-',double(y0),double(v0),'go')
hold off
%
pause;
subplot(1,2,1);
plot(ysol(:,1),ysol(:,2),'r-',double(y0),double(v0),'go')
subplot(1,2,2);
plot(t,ysol(:,1),'b-',t,ysol(:,2),'r-')
legend('y','v')
xlabel 't'
