function fig229
clc
clf
disp('Page 164 Fig2.29-31: Linear System')
disp('dx/dt = 3*x+y')
disp('dy/dt = x-y')
disp('To be shown: vector fields (unnormalized)')
pause;

% Fig. 2.29-31
% linear system
xmin = -3;
xmax = 3;
ymin = -3;
ymax = 3;
dx = 0.3;
dy = 0.3;
[x,y] = meshgrid(xmin:dx:xmax,ymin:dy:ymax);
u = 3*x+y;
v = x-y;
leng = sqrt(u.^2+v.^2);
%
clf
% unnormalized version
quiver(x,y,u,v,0.5)
axis equal tight
xlabel 'x'
ylabel 'y'
title 'Vector Field for a model linear system'
%
disp('To be shown: vector fields (normalized)')
pause
% normalize version
quiver(x,y,u./leng,v./leng,0.4)
hold on
axis equal tight
xlabel 'x'
ylabel 'y'
title 'Vector Field for a model linear system'
%
clc
disp('To be shown: equilibrium points')
pause
% solve for the equilibrium points
[x0,y0] = solve('3*x+y=0','x-y=0');
[x0 y0]
plot(double(x0),double(y0),'go')
xlabel 'x'
ylabel 'y'
%
disp('To be shown: solution curves with initial condition (x,y)(0)=(-0.7082, 3)')
pause;

% solve for the solution curve
fode = inline('[3*y(1)+y(2); y(1)-y(2)]', 't', 'y');
y10 = -0.7082;
y20 = 3;
t0 = 0;
tmax = 4;
[t,ysol] = ode45(fode,[t0 tmax], [y10 y20]);
% solution curve in phase space
plot(ysol(:,1),ysol(:,2),'r-',double(x0),double(y0),'go')
hold off
%
pause;

clf
subplot(1,2,1);
quiver(x,y,u./leng,v./leng,0.4)
hold on
plot(double(x0),double(y0),'go')
subplot(1,2,1);
plot(ysol(:,1),ysol(:,2),'r-')
subplot(1,2,2);
plot(t,ysol(:,1),'b-',t,ysol(:,2),'r-')
legend('x','y')
xlabel 't'
