% here's a sample of how I used matlab to make those figures.  I 
% created this file by first typing "diary for_class".


% this is one of the fixed points.  I found it using maple.
>> a = .1030367663 ; b = .01045779210 ;

% dpdt refers to a file called "dpdt.m" which contains the ODE I'm
% solving.  In your home directory you need a directory called "matlab"
% and in that directory you need the file dpdt.m If your deptartment
% machine doesn't have matlab, come see me.  We have it on the math dept
% pcs and I'll show you how to use it there.

% numerically solve the equation with initial data that is near to (a,b):
>> [t,y] = ode45('dpdt',[0,2000],[a+.1;b]);
% make a plot of x(t) versus t:
>> subplot(3,1,1)
>> plot(t,y(:,1));
>> title('perturbation of steady state (.1030367663, .01045779210)');
>> axis([0,50,-.5,.5])
>> axis([0,50,-.5,.5])
>> subplot(3,1,2)
>> plot(t,y(:,1));
>> axis([0,300,-.5,.5])
>> subplot(3,1,3)
>> plot(t,y(:,1));
% save the figure, to be printed with lpr or ghostview:
>> print -dps fig1.ps

% now I want to make a figure like figure 7.3  So first I go to dpdt.m
% and change k to 0.
>> figure(2)
% figure(2) pops open a new window.  To get to the original window,
% type "figure(1)"
>> a =  -.6734177811; b = 0;
>> [t,y] = ode45('dpdt',[0,500],[a+.1;b]);
>> plot(y(:,1),y(:,2),'-');
>> hold on
% hold on makes it so that I can have more than one orbit per
% figure.  To stop holding on, type "hold off", or "clf" to clear
% the figure-screen altogether.  
>> [t,y] = ode45('dpdt',[0,500],[a+.2;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,500],[a+.3;b]);
>> plot(y(:,1),y(:,2),'-');
% here I noticed that my orbit didn't meet itself so I ran again with
% a longer time interval:
>> [t,y] = ode45('dpdt',[0,700],[a+.3;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,700],[a+.4;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,700],[a+.5;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,700],[a+.6;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,700],[a+.7;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,700],[a+.05;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,700],[a+.005;b]);
>> plot(y(:,1),y(:,2),'-');
% here I decided that I didn't like the curve I just plotted.  So I
% painted over it with itself in black.  (My figure is on a black
% background.)  If my figure was on a white background then I'd do 'w-'
% in the options.  To play around, try options '-', '--', '-.', and ':'.
% These all have color versions 'b-', 'r-', 'y-', 'w-', etc.
>> plot(y(:,1),y(:,2),'k-');
>> [t,y] = ode45('dpdt',[0,700],[a+.02;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,300],[a+1.5;b]);
>> plot(y(:,1),y(:,2),'-');
>> [t,y] = ode45('dpdt',[0,300],[a+1.4;b]);
>> plot(y(:,1),y(:,2),'-');
>> axis([-1.2,1.2,-1.2,1.2])
>> print -dps fig2.ps
>> diary off
>> quit
 4691068 flops.

% note that the above was a little sloppy --- I was always perturbing
% off of the left-most fixed point.  I really should have perturbed off
% of all fixed points and I should have found the eigendirections for
% any saddle points so that I could put initial data on those
% eigendirections to try and capture the unstable and stable manifolds.