Page 40
% make up an m-file, ode243.m, as follows
  % function yprim=ode243(t,y)
  % yprim = y - (y.^2./(2+sin(t)));
% now for the solution with initial value=1
tspan=[0 5];
[t1,y1]=ode23('ode243',tspan,1);
% and for initial value=3
[t3,y3]=ode23('ode243',tspan,3);
% and for initial value=5
[t5,y5]=ode23('ode243',tspan,5);
% plot them all
plot(t1,y1,t3,y3,t5,y5);
 
          Page 43
% Make up an m-file, fig243.m, as follows
  % function Yprime=fig243(t,x)
  % Yprime = [-x(1) - x(2); x(1) - x(2)];
% for the sol. w/ init. value g=1 and i=0
[t,Y]=ode23('fig243',[0 5],[1;0]);
%  semi-colon for column vector
plot(t,Y) % plot both columns of Y vs t

          Page 44
[t,Y4]=ode23('fig243',[0 5],[4;0]);
[t,Y3]=ode23('fig243',[0 5],[3;0]);
[t,Y2]=ode23('fig243',[0 5],[2;0]);
[t,Y1]=ode23('fig243',[0 5],[1;0]);
plot(Y4(:,1),Y4(:,2))
% plot 1st component of Y4 against the 2nd
hold on
plot(Y3(:,1),Y3(:,2)) %ditto for Y3
plot(Y2(:,1),Y2(:,2)) %ditto for Y2
plot(Y1(:,1),Y1(:,2)) %ditto for Y1




          Exercises/Experiments

1.
% to deal with a 2nd order differential
% equation it has to be made into a vector
% valued 1st order DE as follows: the first
% component Y1 is y and the second Y2 is dy/dt.
% Then d2y/dt2+6y=0 becomes the vector system:
%   dY1/dt=Y2; dY2/dt = -6Y1;
% so make an m-file, exer241a.m, as follows
  %  function Yprime=exer241a(t,Y);
  %  Yprime = [Y((2); -6*Y(1)];
[t,Y]=ode23('exer241a',[0 4],[1; 0]);
plot(t,Y(:,1))
%%%
% DE (b) converts to 1st order vector system
%  dY1/dt=Y2; dY2/dt=6*Y1;
% DE (c) converts to 1st order vector system
%  dY1/dt=Y2; dY2/dt=-6*Y1-2*Y2;
% DE (d) converts to 1st order vector system
%  dY1/dt=Y2; dY2/dt=-6*Y1+2*Y2;
% We leave it to the reader to obtain the
% numerical solutions and plot as above.

2a.
% make an m-file, exer242.m, with
%  function Yprime=exer242(t,Y);
%  Yprime=[Y(2); -Y(1)/5+cos(t)];
% then solve and plot with
[t,Y]=ode23('exer242',[0 4*pi],[0;0]);
plot(t,Y(:,1))

2b.
% No built-in direction field in matlab, see
% DFIELD from http://math.rice.edu/~dfield

3a.
% make an m-file, exer243a.m, with
% function yprime=exer243a(t,y);
% yprime=-y.*(1-y);
p=[1 -1 0]; % coef.s of p(y)=-y(1-y)
roots(p)
[t,y]=ode23('exer243a',[0 5],-1);
plot(t,y); hold on
[t,y]=ode23('exer243a',[0 5],-1/2);
plot(t,y)
[t,y]=ode23('exer243a',[0 5],1/2);
plot(t,y)

3b.
% contents of m-file, exer243b.m:
%   function Yprime=exer243b(t,Y);
%   Yprime=[4*Y(1)-Y(1).*Y(1)-Y(1).*Y(2);
%     5*Y(2)-2*Y(2).^2-Y(1).*Y(2)];
[t,Y]=ode23('exer243b',[0 4],[1;1]);
hold off; plot3(t,Y(:,1),Y(:,2))
grid
xlabel('x axis'); ylabel('y axis');
zlabel('z axis'); hold on
[t,Y]=ode23('exer243b',[0 4],[1;4]);
plot3(t,Y(:,1),Y(:,2))
[t,Y]=ode23('exer243b',[0 4],[4;1]);
plot3(t,Y(:,1),Y(:,2))
[t,Y]=ode23('exer243b',[0 4],[4;4]);
plot3(t,Y(:,1),Y(:,2))
view(30,30)   % 30 deg CCW from neg. y-axis, 30 deg elevation
view(-100,30) % 100 deg CW from neg. y-axis, 30 deg elev.

4.
A=[-1, -1; 1, -1]
t=2; At=A*t
expm(At)
t=5; At=A*t
expm(At)
expm(At)*[1;0]
% exp(At)*C where C is a 2x1 column vector