> restart:with(linalg):with(DEtools):with(plots):with(plottools):
# 
# The calculations below are for problem 2 (given eigenvalues and
# eigenvectors construct direction field and general solution).
> d:=matrix(2,2,[-1,0,0,-2]);
> P:=matrix(2,2,[1,-2,1,3]);
> A:=evalm(P&*d&*inverse(P));
> eq:=diff(x(t),t)=-7/5*x(t)+2/5*y(t),diff(y(t),t)=3/5*x(t)-8/5*y(t);
> DEplot([eq],[x,y],t=0..10,x=-5..5,y=-5..5);
# 
# The calculations below are for problem 5 (tank mixing problem).
# 
> eq:=diff(x(t),t)=(-1/2)*x(t)+(1/20)*y(t)+800,diff(y(t),t)=(1/2)*x(t)-(
> 1/4)*y(t),diff(z(t),t)=(1/5)*y(t);
> init:=x(0)=300,y(0)=800,z(0)=0;
> dsolve({eq,init},{x(t),y(t),z(t)});
> assign(");
> subs(t=3,z(t));evalf(");
# 
# The next series of calculations are for problem 6 (body falling in
# water).
# 
> eq6:=m*diff(v(t),t)+10*v(t)=(39/40)*m*g;
> m:=100;g:=981/100;init6:=v(0)=0;
> sol6:=rhs(dsolve({eq6,init6},v(t)));
> solve(sol6-70,t);evalf(");
> 
# 
# The next calculations are for problem 7 (bifurcation diagram).
# 
> eq7:=y^4-7*y^3+14*y^2-8*y-alpha;
> tp1:=textplot({[-4.5,3.3,`(-6.914,3.326)`,align=right],[3,1.53,`(0.94,
> 1.531)`,align=right],[-3.5,0.393,`(-1.383,0.393)`,align=right]}):
> tp2:=textplot({[4,3.8,`sources`],[5,-0.5,`sinks`],[-1.3,1.1,`sources`]
> ,[-4,2.7,`sinks`]}):
> bifur:=implicitplot(eq7,alpha=-8..6,y=-4..8,numpoints=5000,title=`bifu
> rcation diagram`):
> display({tp1,tp2,bifur});
> df:=diff(eq7,y);y_eqpoints:=evalf(solve(df,y));
# 
# 
# Maple has clearly made some numerical roundoff error. The roots we
# have just found should be real numbers. Their imaginary parts are
# zero. So when calculating the corresponding values for alpha, I will
# only use the real parts of these roots.
# 
> i:='i':for i from 1 to 3 do
> solve(subs(y=Re(y_eqpoints[i]),eq7),alpha); 
> od;
> g:=subs(alpha=0,eq7);
> plot(g,y=-2..5,z=-10..10,title=`f(y,0)`);
# 
# The following commands are used to graph the phase line for the
# differential equation when alpha is zero.
# 
> a1:=arrow([0,-2],[0,1.80],.01,.2,.1):a2:=
> arrow([0,5],[0,-4.65],.01,.2,.05):a5:=arrow([0,-2],[0,3.95],.01,.2,.1)
> :a6:=
> arrow([0,5],[0,-2.35],.01,.2,.1):a7:=arrow([0,-2],[0,6.5],.01,.2,.05):
> tplot:=textplot([[1.5,0,`y = 0, sink`],[1.8,1,`y = 1,
> source`,align=right],[1.5,2,`y = 2, sink`,align=right],[1.8,4,`y = 4,
> source`,align=right]]):eqplot:=plot({[[-.1,0],[.1,0]],[[-.1,1],[.1,1]]
> ,[[-.1,2],[.1,2]],[[-.1,4],[.1,4]]},color=black):
> display([eqplot,a1,a2,a5,a6,a7,tplot],scaling=constrained,axes=none,ti
> tle=`Phase line for f(y,0)`);
# 
# The next calculations are for problem 8 (non homogeneous system of two
# equations). 
# 
# 
> eq8:=diff(x(t),t)=-x(t)+2*y(t)-5,diff(y(t),t)=3*x(t)-y(t)-5;
> dfplot:=DEplot([eq8],[x,y],t=0..10,x=-10..10,y=-10..10,color=purple):
> nullclines:=plot([(x+5)/2,3*x-5],x=-10..10,xtickmarks=0,ytickmarks=0,s
> caling=constrained,color=[blue,blue]):
> b1:=arrow([3,4],[4,-2*sqrt(6)],.1,.4,.2,color=red):b2:=arrow([3,4],[4,
> 2*sqrt(6)],.1,.4,.2,color=red):
> display({dfplot,nullclines,b1,b2},title=`Nullclines and Direction
> Field`);
> A8:=matrix(2,2,[-1,2,3,-1]);
> eigenvects(A8);
> evalf(4/3-sqrt(6)/2);
> eq8a:=diff(u(t),t)=-u(t)+2*v(t),diff(v(t),t)=3*u(t)-v(t);
> dsolve({eq8a},{u(t),v(t)});
>