CHAPTER 12
ORDINARY DIFFERENTIAL EQUATIONS

Differential equations provide a bridge between abstract mathmatics, particularly linear algebra and calculus, and their application in dynamic economic models.  The basic problem studied in this chapter is one of finding functions satisfiying a given relationship between a variable and its derivative.  More specifically, the problem is to find function y(t) such that y' (t) = ay(t), for every t in some interval.
Ordinary differential equations contain only ordinary derivatives of order one or higher.  The order of a differential equation is the order of the highest-order derivative that it contains.  The degree of a differential equation is the highest power to which any of the derivatives in the equation is raised.  For example, the differential equation \(y\^\(\^''\)\)(t) + 4\(y\^'\)(t) + y (t) = 0 is of second order, third degree.
A differential equation where the variable and its derivatives are both of the first degree is called a linear differential equation.  For example, the differential equation y' (t) + 2y(t) = 0 is a first order, linear differential equation.

TOOL BOX
Mathematica Commands			Function
DSolve[y' [t] == ay[t] + b, y[t], t]			Solves the differential equation for y(t).
DSolve[{y' [t] == ay[t] + b, y[0] ==yo},  y[t], t]		Solves the initial-value differential equation for y(t).
sol = DSolve[y' [t] == ay[t] + b, y[t], t]			Solves the problem qnd sets the solution equal to sol.
y[t_] = y[t]/.sol[[1]]					uses the replacement rule " /. "to define y[t] as the solu							tion.
NDSolve[{y' [t] == ay[t] + b, y[0]==yo}, y, {t, 0, 10}]  	Finds numerical approximations to the solution.
Plot[Evaluate[y[t]/.%], {t, 0, 10}]			Plots the values of the IntrepolatingFunction.

I.   First order, Homogenous Linear Differential Equations
In many cases, it is observed that the rate of change of some economic variable is proportional to the value of the variable.  For example, the annual investment (the capital growth) of a business may be proportional to its stock of capital, or the aggregate investment in an economy may be proportional to the stock of capital in the economy.  That is,  I(t) = K'(t) = aK(t), where I is the investment and K is the stock of capital.  This, is the simplest form of the differential equations, which is the first-order homogenous linear differential equation.  The problem here is to find the function that gives rise to such a relationship between the variable and its differential.  

In general, the problem of first order homogenous linear differential equations is that of finding the funcition y(t) such that it satisfies, y'(t) = ay(t).  The theory or experiment may supply us with the form of the equation or even the value of the "a".  Any differentiable function that satisfies the above equation is said to be the solution of the differential equation.  To find such a function, divide both sides of the equation y'(t) = ay(t) by y(t) and integrate, ∫((y' (t))/(y (t)))  dt = ∫ a d t.  This will yield, lny(t) + c1 = at + c2 or lny(t) = at + c, where c = c1 + c2.  By transforming this solution into its exponential form we can write, y(t) = (e^c)(e^at).  Set the (e^c) = A, so that the solution will be y(t) = A(e^at).  This solution provides a general solution to the first order, homogenous linear differential equation y'(t) = ay(t).  

In practice, as well as the differential equation, the value of y(t) at some specified time, t = to may be available.  In this case, the set of information, y'(t) = ay(t) and y(to) = yo, is known as the initial-value problem.  Here, the problem is to find a differentiable function y(t), such that y'(t) = ay(t) for all t and y(to) = yo.  Given the general soultion y(t) = A(e^at), y(to) = A(e^ato) =yo or A = yo(e^(-ato)).  Therefore, the initial-value problem y'(t) = ay(t), y(to) = yo has the unique soultion, y(t) = yo(e^(a (t - to))).

In Mathematica, the two commands used to solve the differential equations are "DSolve" and "NDSolve".  The command DSolve solves a differential equation for the function y(t).  The command NDSolve finds numerical approximations to the solution of the differential equation as InterpolatingFunction.  The InterpolatingFunction, then stores a table of values for y(ti) and interpolates this table to find an approximation to y(t) when a particular value of t is given.  The command NDSolve is more appropriate for solving nonlinear differential equations.
 
To solve the initial-value, first order linear differential equation y'(t) = ay(t), y(to) = yo using Mathematica, 
Clear[y]
DSolve[{y'[t]==a y[t], y[to]==yo}, y[t], t]
{{Null (y[t] -> (e^(a t - a to)) yo)}}

Example: Solve y'(t) + 4y(t) = 0, for y(o) = 1.
Solution:  The general soulution of the problems is, y(t) =A(e^(-4 t)).  To find the solution to the initial-value problem, set t = 0 to find A = 1.  Therefore, the soulution of the initial-value problem is y(t) = (e^(-4 t)).
To solve the same problem using Mathematica,
Clear[y]
DSolve[{y'[t]+4y[t]==0, y[0]==1}, y[t], t]
{{Null (y[t] -> (e^(-4 t)))}}

Example:  Using Mathematica, solve y'(t) - .4y(t) = 0, for y(o) = 2.  Find the value of y at t = 5 and plot the time path of y.
Clear[y, sol]
sol = DSolve[{y'[t]-.4y[t]==0, y[0]==2},  y[t], t]
{{y[t] -> 2.` (e^(0.4` t))}}
sol/.t->5
{{y[5]->14.7781}}
Once the problem is solved, use the replacement rule " /. "to define the function y[t] and graph the time path of y.
y[t_]=y[t]/.sol[[1]];
Plot[y[t], {t, 0, 10}]


⁃Graphics⁃

Do More Practice:
1- Change the coefficient of y(t) from -.4 to .4.  Keep the other assumptions the same.  Re-run the program and compare the results.
2- Change the coefficient of y(t) from -.4 to -4.  Keep the other assumptions the same.  Re-run the program and compare the results.
3- Change the coefficient of y(t) from -.4 to  4.  Keep the other assumptions the same.  Re-run the program and compare the results.

Note that, DSolve solves the differential equation for the function y[t].  To find the value of y at t = 5, evaluate the solution as t ->5.   The NDSolve command finds numerical approximations to the solution of the differential equation as InterpolatingFunctionsolution.
Clear[y, sol2]
sol2 = NDSolve[{y'[t] -.4y[t] ==0, y[0]==2}, y, {t, 0, 10}] 
{{y->InterpolatingFunction[{{0.,10.}},<>]}}
NDSolve represents numerical approximations to the solution of the differential equation.  To find a particular value of y at time t, use tthe InterpolatingFunction to estimate y[t].  For example to find the value of y at t = 5, 
y[5]/.sol2
{14.7782}
and to plot the values y[t] for the range t = 0 to t = 10, type
Plot[Evaluate[y[t]/.sol2], {t, 0, 10}]

⁃Graphics⁃

 Do More Practice:  
1- Change the coefficient of y(t) form -.4 to .4.  Solve the problem and observe changes in the solution and the plot.
2- Change the coefficient of y(t) form -.4 to  4.  Solve the problem and observe changes in the solution and the plot.

Example: Application to Business
The net capital investment of a firm is .2 proportion of the capital stock.  Find the time path of investment.  When will the capital stock double?
Solution:  The problem can be written as I(t) = ((dK (t))/dt)= K'(t) = .2k(t).  The solution of this linear homogenous differential equation is K(t) = Ko(e^( .2  t)).  To find the time period after which the capital stock will double set, 2Ko = Ko(e^( .2  t)) or 2 = (e^( .2  t)).  This will yield t = 5ln2 = 3.47.  In Mathematica,
Clear[k]
DSolve[{k'[t] == .2k[t], k[0]==ko}, k[t], t]
{{k[t] -> (e^(0.2` t)) ko}}
Solve[2ko==E^(.20t) ko, t]
{{t->3.46574}}

Do More Practice:
1-  Change the growth rate to .4.  Re-run the program.  How long does it take for capital stock to double at this rate?
2- How long does it take for capital stock to tripple or quadruple at the proportionality rate of .2?

Example: Application in Macroeconomics
Domar's growth model, in a simple form, is a set of equtions as, S = αY, I = β(dy/dt), and S = I.  The first equation expresses actual saving (S) as a constant proportion α of the economy's total output (Y); the second equation assumes investment ( I ) as constant proportion β of the change in total output over time; and the third equation expresses actual saving equal to the economy's desired investment.  Given this information, find the time path of output.
Solution:  Solving the three equations for Y yields, Y = (β/α) (dy/dt), or (dy/dt) = (α/β)y.  The general solution of the problem is Y = A(e^(α(t/β))).  If the level of the output is known for an initail period, say t = to, we can estimate the time path of output for the subsequent periods indefinitely (See,  Evsey D. Domar, "Capital Expansion, Rate of Growth, and Employment, " Econometrica, April, 1946, pp. 137-147; reprinted in Domar, Essays in the Theory of Economic Growth, Oxford University Press, Fair Lawn, N.J., 1957, pp. 70-82.)

Example: Application to Macroeconomics
Suppose that the actual saving of an economy, S(t), is proportional to its total income, Y(t); and the desired investment, I(t), is a constant proportion 4 of the change in income.  Given that marginal propensity to consume (MPC) is .75, and Y(0) = 3000 billion, find the time path of income.  Find y(6).
Solution:  Since MPC = .75, MPS (marginal propensity to save) is .25 and S(t) = .25Y(t).  Also I(t) = 4Y'(t) and Y(0) = 3000.  Using Domar's model above, y'(t) =  .(25/4)y(t).  Using Mathematica, the problem can be solved as
Clear[y, rat]
rat=DSolve[{y'[t]==(.25/4)y[t], y[0]==3000}, y[t], t]
y[t_]=y[t]/.rat[[1]];
Plot[y[t], {t, 0, 10}]
{{y[t] -> 3000.` (e^(0.0625` t))}}

⁃Graphics⁃
y[6]
4364.97

Do More Practice:
1-  Change the MPC to .90 and re-run the program.  Compare the results to those in the text.
2- Suppose that the desired investment, I(t), is a constant proportion .5 of the change in income.  Keep other assumptions the same as in the example.  Re-run the program and compare the results to those in the text.

Problems:
1.  Classify each of the following differential equations as to order and degree.
	a)	y' (t) - y(t) = 4
	b)	ty' (t) + (x^2)(t) - y(t) = 2
	c)	y" (t) - 4y' (t) + y(t) = 0
	d)	\(\(x'\^3\) \((t)\)\) - y' (t) = 1
	e)	(t^3) y" ' (t) + \(\(x'\)'\^2\) = 27
2.   Solve the following differential equations.
	a)	y' (t) + 4y(t) = 0
	b)	y' (t) - 2y(t) = 0
	c)	3y' (t) - 2(e^2)y(t) = 0
	d)	y' (t) - 5y(t) = 0,		y(0) = 10
	e)	3y' (t) - 4y(t) = 0,		y(2) = 4
	f)	y' (t) - 5y(t) = 0,		y(10) = 0		
3.   Suppose that the GDP of a country, y, is growing at a rate proportional to the size of capital stock, k, (dy/dK) = .3k.  Find y as a function of k.		
4.   The population of a coutnry, p(t), is growing at a rate of 1.5 percent a year.  Find the time path of the population growth, given that, population at time t = 0, is 200 million.  How long will it take for the population to double in size?
5.   Suppose that the actual saving of an economy, S(t), is proportional to its total income, Y(t); and the desired investment, I(t), is a constant proportion 3 of the change in income.  Given that marginal propensity to consume is .80, and Y(0) = 200 billion, find the time path of income and y(5).

II.  First Order, Nonhomogenous Linear Differential Equations
A problem of the form y' (t) = ay(t) + f(t), where f(t) is any continuous function of t, is called first order, nonhomogeneous linear differential equation.   Associated with any nonhomogenous equation is the complementary homogenous equation, y'c (t) = ay(t), obtained by setting f(t) = 0.  
		
A solution of the complementary equation is yc(t) = (e^at), the generalization of which provides a solution to the nonhomogeneous equation.  To find the particular solution, yp, search for a function u(t) such that, y(t) = u(t)(e^at).  Diifferentiating this equation yields, y'(t) = au(t)(e^at) + u' (t)(e^at).  Since the first term in the right hand side of this equation, au(t)(e^at), is equal to ay(t), (by the the assumption that y(t) = u(t)(e^at)), the second term, u' (t)(e^at) should be equal to f(t).  Therefore, we can write f(t) = u' (t)(e^at), and we select u(t) such that u'(t) = (e^(-at))f(t).   Having found u(t) by integrating this equation, we obtain the particular solution yp(t) = u(t)(e^at).   However, this solution is only one of many particular solutions to nonhomogeneous equation.  To find the general solution, write y(t) = kyc(t) + yp(t), where k is any constant.  This equation provides a general solution to the nonhomogeneous equation y' (t) = ay(t) + f(t), where yc(t) = (e^at) , yp(t) = u(t)(e^at), and u(t) is obtained by integrating u'(t) = ∫u'(t) dt = ∫ f(t) (e^(-at))dt.  

In the special case where f(t) is a constant and the nonhomogeneous equation is in the form, y' (t) = ay(t) + c, the general solution is, y(t) = k(e^at) - (c/a)(e^at) + (c/a).
	
Example: Solve the nonhomogeneous linear differential equation y' (t) = 2y(t) + (e^(2  t)).
Solution:  The complementary homogenous equation y'(t) = 2y(t) has a solution yc(t) = (e^(2  t)).  Select u(t) such that u' (t) = (e^(-2 t))f(t) = (e^(-2 t)) (e^(2  t)) = 1.  To find u(t), integrate u'(t) in the interval 0 to t to get, u(t) =(∫_0^t) 1 dt = t. Therefore, yp(t) = u(t)(e^at) = t(e^(2  t)), and the general solution to the problem is, y(t) = kyc(t) + yp(t) = k(e^(2  t)) + t(e^(2  t)). Using  Mathematica,
Clear[y]
DSolve[y'[t]==2y[t]+E^(2t), y[t], t]
\({{y[
              t] -> e\^\(2\ t\)\ t + e\^\(2\ t\)\ C[1]}}\)

Do More Practice:
1-  Change the nonhomogeneous part of the problem from (e^(2  t)) to (e^(-2 t)) and re-run the program.  Compare the results to those in the text.
2-  Change the nonhomogeneous part of the problem from (e^(2  t)) to 7 and re-run the program.  Compare the results to those in the text.

Example:  Solve y'(t) = 2y(t) + (e^(2  t)) + 6.
Solution:  The solution follows the solution of the preceding example.  Here, u' (t) = (e^(-2 t))f(t) = (e^(-2 t))((e^(2  t)) + 6) = 1 + 6(e^(-2 t)), and u(t) = (∫_0^t)(1  +  6  (e^(-2 t)) )dt = t - 3(e^(-2 t)) + 3.  The particular solutions is, yp = u(t)(e^at) = (t - 3(e^(-2 t)) + 3)(e^(2  t)) = t(e^(2  t)) + 3(e^(2  t)) - 3, and the general solution of the problem is, y(t) = kyc + yp = k(e^(2  t)) + t(e^(2  t)) + 3(e^(2  t)) - 3 = (k + 3)(e^(2  t)) +  t(e^(2  t)) -3.
Clear[y]
DSolve[y'[t]==2y[t]+E^(2t)+6, y[t], t]
\({{y[
              t] -> \(-3\) + e\^\(2\ t\)\ t + e\^\(2\ t\)\ C[1]}}\)

Example: Application to Markets
Given the following demand and supply functions,  Qd = a - bP, Qs = c + dP, b & d > 0, 
a- Find the equilibrium price and the equilibrium quantity.
b- Suppose that the market price P differs from the equilibrium price Pe and price changes in the market are proportional to excess demand or excess supply in the market.  Find the time path of P.
Solution:  The equilibrium price and quantity are Pe = ((a - c)/(b + d)), and Qe = ((ad + bc)/(b + d)),  respectively.  Assume that the rate of change of price with respect to time is v proportional of excess demand in market, (dP/dt) = v(Qd - Qs), where v is the coefficient of adjustment and (dP/dt) = 0 if Qd =Qs.  Substituting for Qd and Qs from the demand and supply equations we get, (dP/dt) = v(a - bp(t) - c - dp(t)) = -v(b + d)p(t) + v(a - c).  This equation is a first order, nonhomogenous linear differential equation the solution of which is P(t) = k(e^(-v (b + d) t)) + ((a - c)/(b + d)).  Since ((a - c)/(b + d)) = Pe,  P(t) =  k(e^(-v (b + d) t))+ Pe.  If the initial value of the P is given, P(o) = Po, then k = Po - Pe and P(t) = (Po - Pe)(e^(-v(b + d) t))+ pe.  
To solve the same problem using Mathematica, 
ClearAll[Qd, Qs, p]
Qd = a - b p[t];  Qs = c + d p[t]; 
DSolve[{p'[t]==v (Qd - Qs), p[0]==po}, p[t], t]
\({{Null\^2\ \((p[
                  t] -> \(a - c + b\ e\^\(t\ \((\(-b\)\ v - d\ v)\)\)\ \((\(\(-a\) + c\)\/\(b + d\) + po)\) + d\ e\^\(t\ \((\(-b\)\ v - d\ v)\)\)\ \((\(\(-a\) + c\)\/\(b + d\) + po)\)\)\/\(b + d\))\)}}\)
Note that since b, d, and v are assumed to be positive, as t ->∞, (e^(-v(b + d) t)) -> 0.  This implies that, as long as demand and supply have the conventional slopes (negative for demand and positive for supply), the price will always converge to its steady-state equilibrium point.  The steady-state or stationary value of P is the value that (dP/dt) = 0.  However, if demand or supply have unconventional slopes, then the behaviour of the price in market over time will depend on whether b+d is positive or negative.  If b + d > 0,  price will be converge to its steady-state value.  If b + d < 0, then the price will diverge from the steady-state value.  

Example:  Application to Markets
The market demand and supply for a product are given as Qd = 10 - .2p and Qs = 4 + .1p, respectively.   Assume that the price adjustment in the market follows the rule (dp/dt) = .3(Qd - Qs).   Given that the market price at time t = 0 is $4, find the solution of the differential equation and plot the time path of the price.    
ClearAll[p, bat]
Qd=10-2p[t];    Qs=4+.1p[t];
bat=DSolve[{p'[t]==.3(Qd - Qs) , p[0]==4}, p[t], t]
\({{p[t] -> 
            0.14285714285714285`\ e\^\(\(-0.63`\)\ t\)\ \((\(\(8.`\)\(\)\) + 
                  20.`\ e\^\(0.63`\ t\))\)}}\)
p[t_]=p[t]/.bat[[1]];
Plot[p[t], {t, 0, 10}]

⁃Graphics⁃

Do More Practice:
1-  Change the speed of adjustment from .3 to 3.  Re-run the program and compare the results.
2- Change the supply to Qs = -4 + 2p and keep the other assumptions the same as in the text.   Re-run the program and compare the results.
3- Assume that demand has an unconventional slope and is expressed as Qd = 10 + .2p and supply is Qs = 14 + .1p.  Re-run the program and compare the results.
4-  Assume that demand has an unconventional slope and is expressed as Qd = 14+ .1p and supply is expressed as Qs = 10 + .1p.  Re-run the program and compare the results.

Problems:
1. Solve the following differential equations.
	a)	y' (t) + 2 y(t) = 6,		y(0) = 10
	b)	y' (t) + 10 y(t) = 15,		y(0) = 0
	c)	2y' (t) + 3y(t) = 12,		y(0) = 1
	d)	y" (t) - 2y(t) = e^t,		y(1) = 2
	e)	y' (t) - y(t) = t,			y(1) = 10
	f)	y' (t) - y(t) = t + 2,		y(1) = 10
		
2.  Suppose that the GDP growth of country follows, (dy/dk)= .3k + 100.  Find the growth path of y as a function of k.	
3.  The market demand and supply for a product are given as, Qd = 20 - 2p and Qs = 8 + p, respectively.
a)  Find the equilibrium output and price.
b)  Suppose that the market price at time t = 0 is $6.  Find the excess supply in the market.
c)  Assume that the price adjustment in the market follows the rule (dp/dt) = .003(Qd - Qs).  Find the solution of the differential equation and plot the time path of the price.
4.  The demand and supply functions for a product, with a market price of p = 12, are given as, Qd = 700 - 45p and Qs = -100 + 35p, respectively.  Assuming that price adjustment in the market follows the rule (dp/dt) = .01(Qd - Qs).  Find the solution and plot the time path of the price adjustment.
		
III. Systems of First Order, Homogenous Linear Differential Equations
Many economic variables are related to each other through a system of several interrelated variables and equations.  For example, a dynamic IS-LM model of the economy can be expressed by the following system of equations

	(dY/dt) = α[D(Y,r) - Y]		
	(dr/dt) = β[L(Y,r) - M]	
	
Where Y is the output or income, r is the interest rate, M is the money supply, and α and β are the coefficients of adjustment.  The first equation, (dY/dt), represents equilibrium in the output market, IS.  It expresses changes in output as α proportion of the gap between aggregate demand, D(Y, r), and aggregate output, Y.  The second equation,  (dr/dt),  represents the equilibrium in the money market, LM.  It expresses changes in the interest rate as a β propostion of the gap between demand for money, L(Y,r), and the money supply, M.

The laws that govern the behavior of variables in this moded and the system can be estimated by a linear differential equation model.  That is, a model which assumes that the rate of growth of each variable is a linear combination of other variables.  The system in the general form can be written as, Y'(t) = AY(t)  where, Y'(t) is a vector with elements (dYi/dt), A is a matrix of constant coefficients, and Y(t) is the vector of  yi(t)s.

The resemblance of Y'(t) = AY(t) to y'(t) = ay(t) suggets that a trial solution of the form y(t) = (e^at) may be applicable to the system of differential equations.  Let us study the possibility by searching for a constant vector u and a number k such that, y(t) = (e^kt)u satisfies Y'(t) = AY(t).  Since y'(t) = k(e^kt)u, for y(t) to be a solution it is necessary and sufficient that k(e^kt)u = (e^kt)Au.  This equation implies that, Au - ku = (A - kI)u = 0.  Therefore, k is the eigenvalue and u is the eigenvector of the matrix A (see, Chapter 5, Special Topics in Matrix Algebra), and the solution to y'(t) = Ay(t) is, y = (e^kt)u.
	
Example:  Solve the following system of linear homogeneous differential equations, y1' (t) = 2y1(t) + y2(t) and y2' (t) = 3y1(t).
Solution:  The system in the matrix form is, (({{y1' (t)}, {y2' (t)}}))= (({{2, 1}, {3, 0}}))(({{y1 (t)}, {y2 (t)}})).  Find the eigenvalues and eignevectors of A = (({{2, 1}, {3, 0}})),
Eigensystem[(({{2, 1}, {3, 0}}))]
{{-1,3},{{-1,3},{1,1}}}
the two eigenvalues are -1 and 3, with the corresponding eigenvectors (({{-1}, {3}})) and (({{1}, {1}})), respectively.  Therefore, the system of equations have two solutions, y1 = (({{y1 (t)}, {y2 (t)}})) = (e^(-t))(({{-1}, {3}})) and y2 = (({{y1 (t)}, {y2 (t)}})) =  (e^(3  t))(({{1}, {1}})).  
Using Mathematica, the problem can be solved as,
DSolve[{y1'[t]==2y1[t]+y2[t], y2'[t]==3y1[t]}, {y1[t], y2[t]}, t]
{{y1[t] -> 
            (1/4) (e^(-t)) (C[1] + 
                  3 (e^(4 t)) C[1] - 
                  C[2] + (e^(4 t)) C[2]), 
          y2[t] -> 
            (1/4) (e^(-t)) (-3 C[1] + 
                  3 (e^(4 t)) C[1] + 
                  3 C[2] + (e^(4 t)) C[2])}}
For the existence of the solutions of the systems of first order linear differential equations, see appendix 1.

Do More Practice:
1-  Change the coefficient of the y1(t) in the first equation from 2 to -2.  Re-run the program and find the eigenvalues.  What is the effect of the change in the coefficient on the eigenvalues and on the solution.
2-  Change the coefficient of the y1(t) in the first equation from 2 to -2 and the coefficient of the y1(t) in the second equation to -3.  Re-run the program and find the eigenvalues.  What is the effect of the changes in the coefficients on the eigenvalues and on the solution.

Problems:
1.  For the  following matrices, find a solution to y'(t) = Ay(t).
a)  A =  (({{1, -3}, {0, 1}}))		b) A = (({{2, 0}, {1, 0}}))
c)  A = (({{1, 2}, {2, 1}}))			d)  A = (({{2, 3}, {1, 0}}))
e)  A = (({{0, 1}, {(5/3), (2/3)}}))		f)  A = (({{1, 4, 1}, {0, 2, 3}, {0, 0, 3}}))	                 
g)  A = (({{1, 0, 1}, {0, 2, 1}, {4, 0, 1}}))		h)  A = (({{1, 1, 0}, {0, 0, 1}, {0, 1, 1}}))
2.  Solve problems a, b, c, and d of the problem 1 for the initial value, y(0) = (({{1}, {2}})).  Plot the time path of  y(t).
3.  Solve problems a and b of  the problem 1 for the inital value,  y(0) = (({{0}, {3}})).  Plot the time path of  y(t).
4.  Solve problems f and g of the problem 1 for the initial value,  y(0) = (({{1}, {2}, {1}})).
5.  Solve the differential equation  y'(t) = (({{1, 1, 0}, {0, 1, 1}, {0, 0, 1}})) y(t) for the initial value  y(0) = (({{1}, {1}, {1}})).		                  
				
IV. Systems of First Order, Nonhomogenous Linear Differential Equations
A system of first order, nonhomogenous linear differential equations can be expressed in general form as,  y'(t) = Ay(t) + f(t) where f(t) is any continous vector funciton of t.  The resemblance of this function to the first order, nonhomogenous linear differential equation suggests that a solution similar to the solution of y'(t) = ay(t) + f(t) may be applicable here.  An investigation of the solution shows that the equation y(t) = M(t)u(t) + M(t)C is the solution of  y'(t) = Ay(t) + f(t) for any choice of C.  In this solution, IM(t) is the fundamental matrix of the complementary homogeneous system and u(t) is a vector such that,  u'(t) = (M^(-1))(t)f(t).  For definition of the fundamental matrix and the solution of the system of first order nonhomogenous differential equations, the interested reader is referred to appendix 2.  Here, using Mathematica, the application of the system of first order nonhomogenous linear differential equations to economic problems is presented by using several mico and macroeconomic problems. 

Example:  Application to Markets
One of the major characteristics of a perfectly competitive industry is that the firms have easy entry or exit to the industry.  Suppose that the firms' entry or exit in a perfectly competitive industry (the change in the number of firms, (dn/dt)) is proportional to the gap between the market price (p) and the minimum average total cost (atc), (dn/dt) = n' = λ(p - atc), λ > 0.  That is, if p - atc > 0, the positive economic profit attracts more entry to the industry and results in lowers.  The opposite occurs when p - atc < 0.  In addition, assume that changes in the market price is proportional to excess demand in market, (dP/dt) = p' = μ(Qd - Qs),  μ > 0, where Qd and Qs represent market demand and market supply, respectively.  λ and μ represent the speed of adjustment for the number of firms and price, respectively.  Given that Qd = 100 - 2p, Qs =  5n, atc = 20, λ = .1 and μ = .2 formulate the price behaviour in the market, solve the resulting system of differential equations, and plot the time path of the price.  The initial values of p and n are given as po = 25 and no = 15.  Note that, in this problem, it is assumed that the firms in the industry are of the uniform size and that Qs depends on the number of firms in the industry.
Solution:  The system of the differential equations can be expressed as p' = .2(100 -2p -5n) and n' = .1(p - 20).  In matrix form the system can be written as (({{p'}, {n'}})) = (({{- .4, -1}, {.1, 0}}))(({{p}, {n}})) + (({{20}, {-2}})).  The steady-state value of p and n are obtained by setting p' = n' = 0.  Using Mathematica,
Clear[Qd,  Qs,  p,  n,  chimp ]
Solve[(({{- .4, -1}, {.1, 0}})).(({{p}, {n}}))+(({{20}, {-2}}))==(({{0}, {0}})), {p,  n}]
{{Null (p->20.),Null (n->12.)}}
the steady-state values of p and n are 20 and 12, respectively.  To solve the system of differential equations using Mathematica,
Qd = 100 - 2p[t]; Qs = 5n[t];
chimp = DSolve[{p'[t]==.2(Qd-Qs), n'[t]==10(p[t]-20), p[0]==25, n[0]==15}, {p[t], n[t]}, t]
\({{p[
              t] -> \((\(\(0.`\)\(\)\) - 
                  0.010040160642570285`\ i)\)\ \((\((\(\(63.118935352238`\)\(\)\) + 
                        249.`\ i)\)\ e\^\(\((\(-0.2`\) - 3.1559467676119`\ i)\)\ t\) - \((\(\(63.118935352238`\)\(\)\) - 
                        249.`\ i)\)\ e\^\(\((\(-0.2`\) + 3.1559467676119`\ i)\)\ t\) + \((\(\(0.`\)\(\)\) + 
                        1992.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\), 
          n[t] -> \((\(\(0.`\)\(\)\) + 
                  0.002008032128514057`\ i)\)\ \((\((\(\(3992.2726610290533`\)\(\)\) - 
                        747.`\ i)\)\ e\^\(\((\(-0.2`\) - 3.1559467676119`\ i)\)\ t\) - \((\(\(3992.2726610290533`\)\(\)\) + 
                        747.`\ i)\)\ e\^\(\((\(-0.2`\) + 3.1559467676119`\ i)\)\ t\) - \((\(\(0.`\)\(\)\) + 
                        5976.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\)}}\)
p[t_]=p[t]/.chimp[[1]];
Plot[p[t], {t, 0, 10}]

Null (⁃Graphics⁃)

Do More Practice:
1- Change the speed of adjustment for the price equation from .2 to .02.  Keep other assumptions the same.  Re-run the program and compare the results.
2-  Change the speed of adjustment for the number of firms equation from .1 to 10.  Keep other assumptions the same.  Re-run the program and compare the results (Note that, when the speed of adjustment for n is way higher than speed of adjustment for n, at some point, there may be an overshooting of the price.  This happens because the price adjustment depends on the entry and exit of the firms as well as the excess demand).

Problems:		
1.  Solve the following systems of differential equations:
	a)  y'(t) = (({{1, 0}, {-1, 2}}))y(t) + (({{1}, {2}}))					
	
	b)  y'(t) = (({{1, 0}, {-1, 2}})) y(t)+ (({{(e^t)}, {2}}))		
	
	c)   y'(t) = (({{1, 0}, {-1, 1}}))y(t)+ (({{t}, {1}}))	
	
	d)  y'(t) = (({{1, 0, 0}, {0, 2, 0}, {0, 0, 1}}))y(t)+ (({{(e^t)}, {t}, {1}}))	
	
	e)  y'(t)	= (({{1, 0, 1}, {0, 1, 0}, {0, 0, 0}}))y(t) + (({{(e^t)}, {t}, {1}}))					 
2.  Solve problems a, b, and c of the problem for the initial value y(0) =(({{1}, {1}})).		
3.  Solve problems d and e of the problem 1 for the initial value y(0) = (({{1}, {1}, {1}})).		
4.  Suppose that the goods market and the money market of an economy are described by the following system of differential equations:  C = .8Y + 100,  I = .1Y - 200R + 50, Y = C + I,  Md = .8Y - 100R,  Ms = 200,  Md = Ms, where C is consumption, Y is income, I is investment, R is interest rate, Md is money demand, and Ms is money supply.  Output (income) and interest rate adjust to disequilibrium in the market with a speed proportional to excess demand or excess supply in the corresponding markets.  The coefficients of adjustment for the goods market and money market are .01 and .001 respectively.  Build a system of differential equations describing the adjustment of output and interest rate in the markets, and analyze the time path of the output and interest rate. 

V.  Linear Differential Equations of Higher order
The general form of the linear differnetial equations of higher order is, (y^((n))) + an-1(y^((n - 1))) +   ... . ((a1y)^((1))) + aoy = f(t), where (y^((n))) is the nth order derivative of y with respect to t; ais,  i = 0,......, n-1 are constant coefficients; and f(t) is a differentiable function of t.

Solutions of the nth order differential equations are closely related to the soulutions of the systems of differential equations, which were studied in the previous section.  In fact, every nth order equation of the form above can be expressed as a system of differential equations by setting the first order variable (y^((1))) equal to a variable, such as z, which is not differentiated.  The following examples show how this transformtion can be done.

Example:  Solve the second order differential equation \(y\*"""\) - 2y' - 3y = (e^t).
Solution:  Set y' = z.  This implies that, y" = z'  = 2z +3y + (e^t).  These two equations give rise to a system of diffrential equations as, y' = z and z'  = 2z +3y + (e^t).  This system, in a matrix form, can be written as y'(t) = Ay(t) + f(t), where  y'(t) = (({{y'}, {z'}})), A = (({{0, 1}, {3, 2}})), y = (({{y}, {z}})), and  f(t) = (({{0}, {(e^t)}})).  That is, (({{y'}, {z'}}))= (({{0, 1}, {3, 2}}))(({{y}, {z}}))+ (({{0}, {(e^t)}})).  Note that, the second element of the vector y is z, which is the first derivative of y.
The solution of this system of differential equations follows the solution of the systems of first order nonhomogenous linear differential equations of the previous section.  That is, using the two equations u'(t) = (M^(-1))(t)f(t) and y(t) = yp(t) + yc(t)C = M(t)u(t) + M(t)C of the previous section, the solution is y(t) = (({{y}, {z}}))= (1/8)(({{(e^(3  t)) - 2  (e^t) + (e^(-t))}, {3  (e^(3  t)) - 2  (e^t) - (e^(-t))}})).  Using this solution, it can be shown that the second term of the vector y(t) , (z), is the derivative of the first term, (y), and that the solution satisfies the equation y" - 2y' - 3y = (e^t).   Since the initial values are not given, the solution is also one of the general solutions to the probelm. 
To solve the same problem using Mathematica,
ClearAll[y]
DSolve[y''[t] -2y'[t] -3y[t]==E^t, y[t], t]
\({{y[
              t] -> \(-\(e\^t\/4\)\) + e\^\(-t\)\ C[1] + e\^\(3\ t\)\ C[2]}}\)

Example:  Write the third-order differential eqation 2(y^('' '))+ \(\(y\^''\)\) -  y = 2t as a system of differential equations.
Solution:  Set y' = z and y" = z' = x.  This yields the following equations:  y' = z,  z' = x, and x' = -(1/2)x + (1/2)y + t.   In matrix form, these three equations can be written as (({{y'}, {z'}, {x'}})) = (({{0, 1, 0}, {0, 0, 1}, {(1/2), 0, -(1/2)}}))(({{y}, {z}, {x}})) + (({{0}, {0}, {t}})).  The problem can be solved usig the method described in the previous section.  To solve the problem using Mathematica,
ClearAll[y]
DSolve[2y'''[t]+y''[t]-y[t]==2t, y[t], t]
\({{y[
              t] -> \(-2\)\ t + e\^\(t\ Root[\(-1\) + #1\^2 + 2\ #1\^3 &, 1]\)\ C[
                  1] + e\^\(t\ Root[\(-1\) + #1\^2 + 2\ #1\^3 &, 2]\)\ C[
                  2] + e\^\(t\ Root[\(-1\) + #1\^2 + 2\ #1\^3 &, 3]\)\ C[3]}}\)

Companion Matrix:   In general, for every linear differential equation of higher order of the from (y^((n))) + an-1(y^((n - 1))) +   ... . ((a1y)^((1))) + aoy = f(t) there is a corresponding or companion system of first-order linear differential equations, the matrix A of which has the following standard form,  A = (({{0, 1, 0, 0, 0, ., .}, {0, 0, 1, 0, 0, ., .}, {0, 0, 0, 1, 0, ., .}, {., ., ., ., ., ., .}, {., ., ., ., ., ., .}, {-ao, -a1, -a2, -a3, ., ., .}})).   This matrix is called the companion matrix of higher-order equations.  Once the companion matrix of the higher-order differential equation is constructed, the solution of the corresponding system follows the instructions of the previous section.

Example:  Write the companion matrix of the following 4th order differential equation 2  (y^((4))) + 3(y^((2))) -4 (y^((1))) - y = 4(e^t).
Solution:  Write the equation in the standard form , by dividing both sides of equation by 2,   (y^((4))) + (3/2)(y^((2))) - 2 (y^((1))) - (1/2)y = 2(e^t).   The companion matrix of this 4th order differential equation is a 4x4 matrix  A = (({{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {(1/2), 2, ((-3)/2), 0}})).

Example:  Solve the initial-value, second order differential equation y" - 2y' + y = (e^t)  for the initial values y(o) = 1 and y'(o) = 1 ( note that, for the second-order differential equation, the initial value of y and the inital value of the first derivative of y should be given.)
Solution:  The companion system of the problem is y'(t) = (({{0, 1}, {-1, 2}}))y(t) + (({{0}, {(e^t)}})), with the initial values y(0) = (({{1}, {1}})).  The solution to the initial-value problem follows the equation  y(t) = yp(t) + Mn(t)yo,  where yp(t) is the particular solution,  yp(t) = M(t)u(t), M(t) is the fundamental matrix, and Mn(t) is the normalized fundamental matrix,  Mn(t) = M(t)(M^(-1))(0).  To find the fundamental matrix, solve for the eigensystem of the A matrix, 
Eigensystem[(({{0, 1}, {-1, 2}}))]
{{1,1},{{1,1},{0,0}}}
the A matrix has one distinct eigenvalue k = 1 and a distinct corresponding eigenvector u = (({{1}, {1}})).  Therefore, one solution is (({{(e^t)}, {(e^t)}})).  To solve the problem,  we need one more distinct eigenvalue and eigenvector.  This can be obtained from the equation (A - k I)v = u (see the section on the Existence of Definite Solutions).  That is, ((({{0, 1}, {-1, 2}}))- (({{1, 0}, {0, 1}})))(({{v1}, {v2}}))= (({{1}, {1}})).  This results
in v = (({{0}, {1}}))and a second solution for the system as (({{(e^t)}, {(e^t)}}))t + (({{0}, {(e^t)}})) = (({{(te^t)}, {(te^t) +  (e^t)}})).  Therefore, the fundamental matrix of the system is M(t) = (({{(e^t), (te^t)}, {(e^t), (te^t) + (e^t)}})).  Since yp(t) = M(t)U(t) and U'(t) = (M^(-1))f(t), find U'(t) as, 
MatrixForm[Inverse[(({{E^t, t E^t}, {E^t, t E^t +  E^t}}))].(({{0}, {E^t}}))]

(({{-t}, {1}}))
therefore,  U'(t) =  (({{-t}, {1}})).  To find U(t), integrate  U'(t) in the interval 0 to t to find  U(t) = (({{-((t^2)/2)}, {t}})). This will be used to find the particular solution  yp(t) as  yp(t) = M(t)U(t). 
MatrixForm[(({{(e^t), t (e^t)}, {(e^t), t (e^t) + (e^t)}})).(({{-((t^2)/2)}, {t}}))]
(({{(((e^t) (t^2))/2)}, {-(1/2) (e^t) (t^2) + t ((e^t) + (e^t) t)}}))
To find the initial value solution, first find the normalized matrix Mn(t) = M(t)(M^(-1))(0),
MatrixForm[Inverse[(({{E^t, t  E^t}, {E^t, t  E^t + E^t}}))]/.{t -> 0}]
(({{1, 0}, {-1, 1}}))
MatrixForm[(({{E^t, t  E^t}, {E^t, t  E^t + E^t}})).(({{1, 0}, {-1, 1}}))]
(({{(e^t) - (e^t) t, \(e\^t\ t\)}, {-(e^t) t, \(e\^t + e\^t\ t\)}}))
 M(t)(M^(-1))(0)= (({{(e^t) - (te^t), (te^t)}, {-(te^e), (e^t) + (te^t)}})).  To find  y(t) = yp(t) + Mn(t)yo, 
MatrixForm[(({{(((e^t) (t^2))/2)}, {-(1/2) (e^t) (t^2) + t ((e^t) + (e^t) t)}}))+(({{(e^t) - t (e^t), t (e^t)}, {-t (e^t), (e^t) + t (e^t)}})).(({{1}, {1}}))]
(({{(e^t) + (((e^t) (t^2))/2)}, {(e^t) - (((e^t) (t^2))/2) + t ((e^t) + (e^t) t)}}))
Therefore, y(t) = (e^t) + (1/2) (e^t) (t^2).  
To solve the same initial-value problem using Mathematica,
ClearAll[y]
DSolve[{y''[t] -2y'[t]+y[t]==E^t, y[0] == 1, y'[0] ==1}, y[t], t]
{{y[t] -> (e^t) (1 + ((t^2)/2))}}

Example:  Application to Markets
The market demand and supply of a good produced and sold in a perfectly competitive market are given as Qdt = 100 - .2Pt and Qst = -10 + .3Pet, respectively.  P represents the market price and Pe represents the expectation of the market price by the producers.  It is assumed that producers form their expectations of market prices by Pet = Pt + .4P't + .2P''t rule, where P't = (dP/dt) and P''t = (((d^2) P)/(dt^2)).   In addition, changes in the market price is proportional to excess demand in market with a proportionality factor of .3.  That is (dP/dt) = .3( Qd - Qs).  Given that the initial prices are given P(0) = 200, P(1) = 210, solve the dynamic market model and plot the time path of Pt.
Solution:  Substitute for Qd and Qs in the price adjustment equation (dP/dt) = .3( Qd - Qs) to get (dP/dt) = .3(100 - .2P +10 -.3Pet).  Susbstitue for Pet  to get .15Pt +1.036P't + .018P''t = 33.  This is a second order, nonhomogeneous differential equation.  The steady-state value of the price can be derived by setting p' = p'' = 0.   The steady-state value of the price is P = 220.  The differential equation can be solved using Mathematica as,
Clear[Qd, Qs, p, koala]
Qd = 100-.2p[t]; Qs = -10+.3(p[t]+.4p'[t]+.2p''[t]);
koala = DSolve[{p'[t]==.3(Qd - Qs), p[0]==200, p[1]==210}, p[t], t]
\({{Null\^2\ \((p[
                  t] -> e\^\(\(-0.14515371895586968`\)\ t\)\ \((\(-11.562172891227148`\) - 
                      8.437827108772852`\ e\^\(\(-57.26524811764381`\)\ t\) + 
                      220.`\ e\^\(0.14515371895586968`\ t\))\))\)}}\)
p[t_]=p[t]/.koala[[1]];
Plot[p[t], {t, 0, 50}]

Null (⁃Graphics⁃)

Do More Practice:
1-  Change the speed of price adjustment to excess demand from .3 to 3.  Keep other assumptions thesame.  Re-run the program and compare the results.
2- Change the speed of price adjustment to excess demand from .3 to .03.  Keep other assumptions thesame.  Re-run the program and compare the results.
3- Change the Pe equation to  Pet = Pt + .6P't + .1P''t.  Keep the other assumptions the same.  Re-run the program and compare the results.	

Example:  Application to Markets
In the section for the First Order, Nonhomogeneous Linear Differential Equations, the market demand and supply for a product were given as, Qd = 10 - .2p and Qs = 4 + .1p, respectively.  It was assumed that the price adjustment in the market was proportional to the excess demand in market.   Now, suppose that the price adjustment in the market depends on the accumulation of the stock of unsold inventories of the product as well as the excess demand.  That is, the excess demand in market inceases the market prices by the rule (dp/dt) = .3(Qd - Qs) and the stock of unsold inventories depresses the market prices by the rule (dp/dt) = -.08(∫_0^t)(Qs - Qd) dt.   Given the initial market prices of Po = $6 and P1 = 5.8, find the solution of the differential equation and plot the time path of the price.    
Solution:  The problem can be formulated as (dp/dt)=  .3(Qd - Qs) - .08(∫_0^t)(Qs - Qd) dt.  Differentiate the function with respect to t to get (((d^2) p)/(dt^2))=  .3((d(Qd  -  Qs))/dt) - .08 (Qs - Qd).  Substitute for Qd and Qs from the demand and supply equations and differentiate with respect to t to get  P'' = -.09P' - .024P + .48, where P'' = (((d^2) p)/(dt^2)) and P' = (dp/dt).   The steady-state value of P is derived by setting P'' = P' = 0.  This yields P = 20.  To find the adjustment path of P over time solve the second order, nonhomogeneous differential equation above and analyze the time path of Pt.  Using Mathematica, 
Clear[p, ox]
ox=DSolve[{p''[t]==-.09p'[t]-.024p[t]+.48, p[0]==6, p[1]==5.8}, p[t], t]
\({{Null\ \((p[
                  t] -> e\^\(\(-0.045`\)\ t\)\ \((20.`\ e\^\(0.045`\ t\) - 14.`\ Cos[0.14823967080373593`\ t] - 
                      6.818938852372753`\ Sin[
                          0.14823967080373593`\ t])\))\)}}\)
p[t_]=p[t]/.ox[[1]];
Plot[p[t], {t, 0, 100}]

Null (⁃Graphics⁃)
It is obvious from the graph that P converges to its steady state equilibrium value with an oscillatory path.  

Do More Practice:
1-  Change the speed of price adjustment to excess demand from .3 to 3.  Keep other assumptions thesame.  Re-run the program and compare the results.
2- Change the speed of price adjustment to excess demand from .3 to .03.   Change the speed of price adjustment to accumulation of inventories from .08 to .3.  Keep other assumptions thesame.  Re-run the program and compare the results.
3- Change the supply equation to Qs = -4 + 2p and keep the other assumptions of the problem the same.   Re-run the program and compare the results.	

Problems:
1.  Write the companion systems of the following differential equations.
	a)	(y^((2))) - (y^((1))) - 2y = 0
 	b)	3(y^((2))) - 2(y^((1))) - 5y = 0
 	c)	(y^((3))) - y = 0
 	d)	(y^((3))) - (y^((2))) + 2(y^((1))) - 3y = (e^t)
 	e)	(y^((3))) + 2(y^((2))) - 5y = (e^t)
 	f)	3(y^((4))) - 6(y^((2))) + 3(y^((1))) + y = 0 
 2.  Solve the following differential equations.
	a)	(y^((2))) - 2(y^((1))) + y = (e^(2  t))
	b)	(y^((2))) - 2(y^((1)))+ y = 0
	c)	3(y^((2))) + 2(y^((1))) - y = 0
	d)	(y^((2))) - 4Y = 3(e^t)
	e)	3(y^((2))) - 2(y^((1))) - 5y = 0 		
3.  Solve equation a of the problem 2 for the initial values y(o) = 1,  y'(o) = 2.
4.  Solve equation b of the problem 2 for the initial values y(o) = 0,  y'(o) = 1.
5.  Solve equation c of the problem 2 for the initial values y(o) = 1, y'(o)= 1.
6.  Solve equation d of the problem 2 for the initial values y(o) = 1, y' (o) = 1.
7.   The market demand and supply for a product are given as Qd = 50 - .2p and Qs = -10 + .3p, respectively.  It is assumed that the price adjustment in the market is proportional to excess demand in market, with a proportionality factor of .4.  As well, it is assumed that the price adjustement in market is proportional to the accumulation of the stock of unsold inventories, with a proportionality factor of .1.   Given the initial prices of Po = 100 and P1 = 105, find the steady-state value of the price, solve the dynamic model of the market for Pt, and plot the time path of Pt.   Analyze the adjustment process.  
 
VI. Differential Equations and Dynamic Stability of Economic Models
Economic models, in general, consist of a set of behavioral relations and a set of equilibrium conditions.  For example, a simple market model consists of two behavioral equations, Qd = a - bp , Qs = c + dp,  and a market-clearing or equilibrium condition, Qd = Qs.  In these equations, Qd and Qs are quantity demanded and quantity supplied respectively, p is the market price, a, b, c, and d are parameters, where b & d > 0.  Solving the system for the equilibrium condition yields pe = ((a - c)/(b + d))and Qe = ((ad + cb)/(b + d)), where Qe and pe are equilibrium quantity and equilibrium price, respectively.
 
The equilibrium values derived above represent only the solution to the system of equations.  However, the solution is a static one and does not say anything about the time path by which the system moves from one equilibrium to another.  The differential equations of this chapter are dynamic models, where time is explicitly included in the models, and the models specify how variables react when they are away from their equilibrium values, or how they move from one equilibrium to another one.  The stability analysis of economic models in this section is the study of how the equilibrium values of the variables react if there is any disturbance in the market.
 
 An equilibrium may be stable, unstable, or neutral.  An equilibrium is called stable if the values of the variables in question return to the equilibrium following a disturbance.  An equilibrium is called unstable if the values of the variables progressively diverge from the equilibrium point as time passes.  An equilibrium is neutral if the values of the variables remain unchanged.
 
 In this section we will study the stability of economic models using first-order and higher-order linear differential equations.
 
 Stability Analysis of the First-Order Linear Homogeneous Differential Models:   An initial-value first order linear homogenous differential equation y'+ ay = 0,  y(0) = yo,  has a single (characteristic) root, -a, and the general solution of the equation is, y(t) = yo(e^(-at))	.  As t -> ∞,  it is clear that the value of y will diverge from yo.   If a is positive, as  t -> ∞,  y will approach zero.  if a is a negative number, as  t -> ∞,  y will grow exponentially away from yo and will approach infinity.  Therefore, the stability of an economic model with a first-order linear differential equation will depend on the sign of the (characteristic) root of the differential equation.
 
 Example:  Solve the initial-value differential equation y' + 2y = 0, y(0) = 1and analyze the stability of the model.
  Solution:  Since a = 2 > 0, as  t -> ∞, the value of y will converge to zero.  To show this graphically, solve the problem and plot the time path of y,
ClearAll[y, fox]
 fox=DSolve[{y'[t]+2y[t]==0, y[0]==1}, y[t], t]
{{y[t] -> (e^(-2 t))}}
y[t_]=y[t]/.fox[[1]]
Plot[y[t], {t, 0, 5}]
(e^(-2 t))

⁃Graphics⁃

Do More Practice:  
Change the value of "a" to -2.  Re-run the program and observe changes in solution and the plot.
 
 The conclusion derived in the stability analysis above is true also for the initial-value, first-order nonhomogeneous differential equations of the form y' + ay = b.   If the value of "a" in the equation y' + ay = b is positive, as t -> ∞,  y will approach zero.  if a is a negative number, as  t -> ∞,  y will grow exponentially away from yo and will approach infinity.  However, if the nonhomogeneous part of the differential equation is a function of t, y' + ay = f(t), then the stability of the model will depend on the form of the f(t) function.   For these modes, the best way to analyze the stability of the model is to plot the time path of y and observe its bebavior over time.
 
 Example:  Solve the initial-value differential equation y' + 2y = 3, y(0) = 1and analyze the stability of the model.
ClearAll [y, dog]
 dog=DSolve[{y'[t]+2y[t]==3, y[0]==1}, y[t], t]
{{y[t] -> 
            (1/2) (e^(-2 t)) (-1 + 
                  3 (e^(2 t)))}}
y[t_]=y[t]/.dog[[1]]
Plot[y[t], {t, 0, 5}]
(1/2) (e^(-2 t)) (-1 + 
            3 (e^(2 t)))

⁃Graphics⁃

 Do More Practice:  
1- Change the value of "a" to -2.  Re-run the program and observe changes in solution and the plot.
2- Change the value of  "b"  to -3.  Re-run the program and observe changes in solution and the plot.

Example:  Solve the initial-value differential equation y' + 2y = 3t, y(0) = 1and analyze the stability of the model.
Clear[y, bee]
bee=DSolve[{y'[t]+2y[t]==3t, y[0]==1}, y[t], t]
{{Null (y[t] -> 
                (1/4) (e^(-2 t)) (7 - 
                      3 (e^(2 t)) + 
                      6 (e^(2 t)) t))}}
y[t_]=y[t]/.bee[[1]]
Plot[y[t], {t, 0, 5}]
(1/4) (e^(-2 t)) (7 - 
            3 (e^(2 t)) + 
            6 (e^(2 t)) t)

⁃Graphics⁃

 Do More Practice:  
1- Change the value of "a" to -2.  Re-run the program and observe changes in the solution and the plot.
2- Change the value of  f(t) to -3t.  Re-run the program and observe changes in the solution and the plot.
3- Change the value of  f(t) to 3(e^t).   Re-run the program and observe changes in the solution and the plot.
4- Change the value of  f(t) to 3(e^(-t)).  Re-run the program and observe changes in the solution and the plot.

 Example:  Application to Markets
 Suppose that in the market model above the market price p ≠ pe.   Generally, price changes in the market are proportional to excess demand in the market, i.e. (dp/dt) = λ(Qd - Qs), λ > 0, and (dp/dt) = 0 only if Qd = Qs.  Here, λ is the coefficient of adjustment or speed of adjustment .  Substituting for Qd and Qs from the demand and supply equations, (dp/dt) = λ(a -c) - λ(b + d)p,  or (dp/dt)+ λ(b + d)p = λ(a - c).   This is a first-order linear differential equation, the solution of which is, p(t) = (p(o) - ((a  -  c)/(b + d)))(e^(-λ (b + d) t)) + ((a - c)/(b + d .  ))Since pe = ((a - c)/(b + d)), we can write p(t) = (p(o) - pe)(e^(-λ (b + d) t)) + ((a - c)/(b + d)).   Therefore, the dynamic stability of the market requires that, λ(b+d) > 0.   Since λ, b, and d are positive by assumption, the stability condition is satisfied as long as demand and supply have the conventional downward sloping and upward sloping shapes, respectively.   However, if demand or supply have unconventional slopes, the stability condition would require that b > -d.  That is, for the model to be stable, the slope of the supply curve should exceed the slope of the demand curve.  If this condition is satisfied,  (e^(-λ (b + d) t)) will approach zero and p(t) will approach pe as t -> ∞.  Otherwise, p(t) will progressively deviate from pe.
 
 Example:  Application to Markets
 Demand and supply for a product are given as Qd = 10 - .2p and Qs = -2 + p, respectively.  A shift in demand has resulted at a market price of p = 8.  Suppose that the price changes in the market adjust to excess demand with an adjustment rate of  (dp/dt)= .04(Qd - Qs).   Find the equilibrium price and equilibrium quantity.  Find the excess demand at p = 8 and the time path of the price adjustment to equilibrium.  Analyze the stability of the market.
  Solution:  To find the equilibrium price and quantitiy set Qd = Qs and solve for p and Q,
ClearAll[p, owl] 
Qd = 10 - .2p[t];  Qs = -2 + p[t];
 owl = Solve[Qd==Qs, p[t]]
{{Null (p[t]->10.)}}
Qd/.owl
{8.}
Therefore, the equilibrium price and quantity are $10 and 8, respectively.  To find the excess demand at p = 8,
(Qd-Qs)/.p[t]->8  
2.4
the excess demand is 2.4 units.  To find the time path of the price adjustment to equilibrium, solve the differential equation  (dp/dt) = .04(Qd - Qs), p(o) = 8, and plot the time path of the price,
ClearAll[p, owl2]
owl2 = DSolve[{p'[t]==.04(Qd-Qs), p[0]==8}, p[t], t] 
{{p[
              t] -> (e^(-0.048` t)) (-2.` + 
                  10.` (e^(0.048` t)))}}
p[t_]=p[t]/.owl2[[1]]
Plot[p[t], {t, 0, 100}]
(e^(-0.048` t)) (-2.` + 
            10.` (e^(0.048` t)))

⁃Graphics⁃
the time path of the price adjustment shows that as t -> ∞, the market price approaches the equilibrium price of p = 10.

 Do More Practice:  
1- Change the speed of adjutment  λ in the (dp/dt) = λ(Qd - Qs) from .04 to a larger or smaller number, (.4 or .004).  Re-run the program and observe the cahnges in the solution and the time path of the price adjustment.   
2- Change the demand and supply equations of the preceding example to Qd = 10 - .8p and Qs = 4 - .2p.   Keep the equation of the adjustment and other assumptions constant.  Solve the resulting differential equation and graph the time path of the price adjutment by re-running the program.
3- Change the demand and supply equations of the preceding example to Qd = 10 - .2p and Qs = 16 - .8p.   Keep the equation of the adjustment and other assumptions as it is.  Re-run the program and compare the solution and the graph to the one in the text.
 
Stability Analysis of Linear Homogeneous Differential Models of Higher Order:   A second-order linear differential equation y" + ay' + by = 0 has a companion system y' = A y =(({{0, 1}, {-b, -a}}))y.  The characteristic roots, or the eigenvalues of the A matrix of the system are a solution to the quadratic equation (k^2) + ak + b = 0 with the roots k1 = (1/2) (-a  +  (Sqrt[(a^2)  - 4  b])) and k2 = (1/2) (-a  -  (Sqrt[(a^2)  - 4  b])).   If k1 and k2 are distinct, the system has two distinct characteristic vectors, u1 and u2, corresponding to the characteristic roots k1 and k2.  Also, the system will have two distinct solutions,  y1 = u1(e^k1t) and y2 = u2(e^k2t).  As well, any linear combination of y1 and y2 such as,  y = c1y1 + c2y2 = c1u1(e^k1t) + c2u2(e^k2t) will also be a solution to the system.  The values of c1 and c2, however, are determined by the initial conditions.  

If (a^2) - 4b of the characteristic roots above is positive, the roots are real and it is clear, from the solution y = c1y1 + c2y2 = c1u1(e^k1t) + c2u2(e^k2t)  that the time path of  y  will approach equilibrium, if k1 and k2 are both negative.  For the characteristic roots k1 and k2 to be negative, the solutions k1 = (1/2) (-a  +  (Sqrt[(a^2)  - 4  b])) and k2 = (1/2) (-a  -  (Sqrt[(a^2)  - 4  b]))  imply that a and b both have to be positive.  Therefore, for the model y'' + ay' + by = 0 to be stable, both a and b have to be positive.
 
If only one of the characteristic roots is positive while the other one is negative, say for example if k1< 0 and k2 > 0, then the term (e^k1t) in the solution equation will decay to zero as  t -> ∞,  while the term (e^k2t) will grow exponentially away from the equilibrium.  In this case, the model has a saddle point.

If (a^2) - 4b of the characteristic roots is negative, the roots are complex conjugates.   If the real parts of the complex conjugate roots are negative, the time path of y will spiral into equilibrium.   The solution and the stability of equations with the complex conjugate roots will be discussed in chapter 14.
 
 In general, for every nth-order linear homogeneous differential equaion  (y^((n))) + an-1(y^((n - 1))) +   ... . ((a1y)^((1))) + aoy = 0  there is a characteristic equation of the form  (k^(n ))+ an-1(k^(n - 1)) + . . . . . . + a1k + ao = 0  the solution of which gives the characteristic roots of the differential equation.  The stability of the nth order differential equations requires that all the ai coefficients of the characteristic equation to be positive.  If one or more of the ai coefficients are negative, then the system possesses positive roots, and it is dynamically unstable.
 
Eigenvalues and the Stability Analysis of the Linear Homogeneous Differential Equations:  The stability test of the linear differential equations can also be done using the eigenvalue test on the A matrix of the system of linear differential equations or on the companion matrix of the companion system of higher order linear differential equations.  A system of linear homogeneous differential equations or a higher order homogeneous linear differential will be stable iff the real parts of the eigenvalues of the A matrix of the system or the real parts of the eigenvalues of the companion matrix are all negative.   If the eigenvalues are real with some positive and some negative, we have a saddle point.  If the eigenvalues are complex conjugates with negative real parts, all paths spiral into the equilibrium (stability of the systems with complex conjugate eigenvalues will be studied in chapter 14).  

The above conditon is equivalent to saying that the necessary and sufficient conditins for the stability of the nth order differential equations are that the trace of the corresponding companion matrix to be negative and the determinant of the corresponding companion matrix to be positive.  In many occasions, the trace and determinant conditions are a more convenient way of analyzing the stability of some higher order linear differential economic models.
 
 Example:   Analyze the stability of the second-order differential equation 	y" + 5y' + 6y = 0, y'(0) = y(0) = 1.
  Solution:   The corresponding characteristic equation of the system is (k^2) + 5k + 6 = 0.   Since all the coefficients of the characteristic equation are positive, the characteristic roots are negative and the system is dynamically stable.  The characteristic roots of the equation are -2 and -3.
  A second way to do the same analysis is to form the companion matrix of the corresponding companion system, A = (({{0, 1}, {-6, -5}})).   The trace of the matrix is tra(A) = -5 and the determinant is det(A) = 6.  Since the trace is negative and the determinant is positive, the necessary and sufficient conditons for the stability of the system are satisfied and the system is dynamically stable.
  To solve for y and to analyze the stability of the model using Mathematica,
ClearAll[y, elk] 
elk=DSolve[{y''[t]+5y'[t]+6y[t]==0, y'[0]==y[0]==1}, y[t], t]
{{y[
              t] -> (e^(-3 t)) (-3 + 
                  4 (e^t))}}
 y[t_]=y[t]/.elk[[1]];
 Plot[y[t], {t, 0, 10}]

⁃Graphics⁃

 Do More Practice:  
1- Change the coefficients form 5 and 6 to 5 and -6 and observe changes in the solution and the plot.
2- Change the coefficients to 3 and 6 and observe changes in the solution and the plot.

Example:  Solve the following system of linear homogeneous differential equations and analyze the stability of the model, Y' = (({{3, 1}, {1, 2}}))y.
Solution:  The eigenvalues of the A matrix are,
Eigenvalues[(({{3, 1}, {2, 1}}))]
{(2 - (Sqrt[3])) Null, (2 + (Sqrt[3])) Null}
the eigenvalues are both real and positive.  Therefore, the system will be an unstable one.
To solve the problem using Mathematica,
ClearAll[y1, y2, hen]
hen=DSolve[{y1'[t]==3y1[t]+y2[t], y2'[t]==y1[t]+2y2[t], y1[0]==y2[0]==1}, {y1[t], y2[t]}, t]
\({{y1[t] -> 
            1\/10\ \((5\ e\^\(1\/2\ \((5 - \@5)\)\ t\) - 
                  3\ \@5\ e\^\(1\/2\ \((5 - \@5)\)\ t\) + 
                  5\ e\^\(1\/2\ \((5 + \@5)\)\ t\) + 
                  3\ \@5\ e\^\(1\/2\ \((5 + \@5)\)\ t\))\), 
          y2[t] -> 
            1\/10\ \((5\ e\^\(1\/2\ \((5 - \@5)\)\ t\) - \@5\ e\^\(1\/2\ \((5 - \@5)\)\ t\) + 
                  5\ e\^\(1\/2\ \((5 + \@5)\)\ t\) + \@5\ e\^\(1\/2\ \((5 + \@5)\)\ t\))\)}}\)
y1[t_]=y1[t]/.hen[[1]]
y2[t_]=y2[t]/.hen[[1]]
(1/10) (5 (e^((1/2) (5 - (Sqrt[5])) t)) - 
            3 (Sqrt[5]) (e^((1/2) (5 - (Sqrt[5])) t)) + 
            5 (e^((1/2) (5 + (Sqrt[5])) t)) + 
            3 (Sqrt[5]) (e^((1/2) (5 + (Sqrt[5])) t)))
\(1\/10\ \((5\ e\^\(1\/2\ \((5 - \@5)\)\ t\) - \@5\ e\^\(1\/2\ \((5 - \@5)\)\ t\) + 
            5\ e\^\(1\/2\ \((5 + \@5)\)\ t\) + \@5\ e\^\(1\/2\ \((5 + \@5)\)\ t\))\)\)
the plot of y1 and y2 over time is shown below,
Plot[y1[t], {t, 0, 5}]
Plot[y2[t], {t, 0, 5}]


(⁃Graphics⁃) (⁃Graphics⁃)
the parametric plot of y1, y2 , and t is shown below,
ParametricPlot3D[{t, y1[t], y2[t]}, {t, 0, 10}, Boxed->False, BoxRatios->{5, 2, 2}, AxesLabel->{t, y1, y2}]


⁃Graphics3D⁃
all three graphs show that the system is unstable.

 Do More Practice:  
1- Change the coefficients of the A matrix to (({{1, 0}, {1, 4}})) and re-run the program.  Compare the changes.
2- Change the coefficients of the A matrix to (({{0, 1}, {3, 4}})) and re-run the program.  Compare the changes.
3- Change the coefficients of the A matrix to (({{1, -1}, {4, 4}})) and re-run the program.  Compare the changes.

Example: Application to IS-LM Model
Suppose that in an IS-LM model of the economy, income(Y) adjusts with an α proportion to the excess demand in the product market and interest rate (r) adjusts with a β proportion to the excess demand for money in the money market.  That is, y' = α(D(y, r) - y) and r' = β(L(y, r) - M), where D(y, r) is aggregate demand, L(y, r) is demand for money, and M is money supply.  Analyze the stability model with respect to the slopes of the IS and LM and the speeds of adjustments in the two markets.
 Solution:  Linearizing the model (totally differentiating the functions) in the neighborhood of the equilibrium point , y = ye and r = re, yields  dy'  = α(Dydy + Drdr - dy) and dr'  = β(Lydy + Lrdr - dM).  In the matrix form the system can be written as (({{dy'}, {dr'}}))= A(({{dy}, {dr}})) =(({{-α (1 - Dy), αDr}, {βLy, βLr}}))(({{dy}, {dr}})), where it is assumed that money supply is decided exogenously by the Fed and dM = 0.
The stability of the system requires that the trace of the A matrix of the system to be negative and its determinant to be positive.  That is tr(A)   = -α(1 - Dy) + βLr < 0 and det(A) = - αβLr(1 - Dy) - αβLyDr > 0.  The trace condition will be satisfied if (1 - Dy) > 0 and Lr < 0.   With the conventional assumptins on the slopes of the IS and LM curves, both of these conditions are satisfied.  The determinant conditon for the stability requires that ((1  -  Dy)/Dr) < - (Ly/Lr).   Here, ((1  -  Dy)/Dr) is the slope of the IS curve and - (Ly/Lr) is the slope of the LM curve.  The stability of the system requires that the slope of the LM curve exceeds the slope of the IS curve (Note:  To find the slope of the IS and LM curves, totally differentiate IS and LM equations at the equilibrium point, where D(Y, r) = y and L(y, r) = M, then find (dr/dy) for each equation.  This will yield (dr/dy) = ((1  -  Dy)/Dr)as the slope ofthe IS curve, and (dr/dy) = - (Ly/Lr) as the slope of the LM curve.	    
	
Example:  Application to IS-LM Models
A simple macro model of an economy is represented by the following equations: c = 500 + .9y,  i = 100 - 3r, y = c + i,  Md = .25y - .8r, and Md = Ms = 400.  The first equation expresses consumption (c) as a function of income (y); the second equation expresses investment ( i ) as a function of interest rate (r); the third equation is the equilibrium condition in the goods market; the fourth equation expresses demand for mone (Md) as a function of income and interest rate; and the last equation is the equilibrium condition in the money market.  Suppose that the speed of the adjustments in the goods market and money market are proportional to the excess demand in the corresponding markets, with a proportionality factor of .4 for the goods market and a proportionality factor of .6 for the money markets.  Given that the initial values of y and r are y(0) = 3000 and r(0) = 10, find the time path of y and r and analyze the stability of the model.  
Solution:  The equilibrium condition in the goods market requires y = c + i = 600 + .9y - 3r and the equilibrium condition in the money market requires .25y - .8 r = 400.  The slope of the IS curve is, (dr/dy)=  .(10/(-3)) = - .033 and the slope of the LM curve is (dr/dy)= \(TraditionalForm\` .25\/ .8\)= .312.   IS and LM have the conventional slopes and the trace condition for the stability of the model is satisfied.  The determinant condition for the stability of the model requires that the slope of the LM curve exceeds the slope of the IS curve.  This condition is also satisfied.   Therefore, the model is a stable model.  
To solve the differential equation using Mathematica,
ClearAll[y,c, i,  r,Md, Ms,  hog]
c = 500+.9y[t]; i = 100 - 3r[t]; Md = .25y[t]-.8r[t]; Ms = 400;
hog=DSolve[{y'[t]==.4(c + i -y[t]), r'[t]==.6(Md-Ms), r[0]==10, y[0]==2000}, {y[t], r[t]}, t]
\({{Null\^2\ \((y[
                  t] -> \((\(\(0.`\)\(\)\) + 
                      0.0036620646720621083`\ i)\)\ \((\((\(\(53344.90836996536`\)\(\)\) + 
                            3290.`\ i)\)\ e\^\(\((\(-0.26`\) - 0.3627671429443411`\ i)\)\ t\) - \((\(\(53344.90836996536`\)\(\)\) - 
                            3290.`\ i)\)\ e\^\(\((\(-0.26`\) + 0.3627671429443411`\ i)\)\ t\) - \((\(\(0.`\)\(\)\) + 
                            552720.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\))\), 
          Null\^2\ \((r[
                  t] -> \((\(\(0.`\)\(\)\) + 
                      0.00018310323360310542`\ i)\)\ \((\((\(\(175706.2656850916`\)\(\)\) + 
                            334593.`\ i)\)\ e\^\(\((\(-0.26`\) - 0.3627671429443411`\ i)\)\ t\) - \((\(\(175706.2656850916`\)\(\)\) - 
                            334593.`\ i)\)\ e\^\(\((\(-0.26`\) + 0.3627671429443411`\ i)\)\ t\) - \((\(\(0.`\)\(\)\) + 
                            723800.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\))\)}}\)
To show the dynamics of the time path to equilibrium, 
y[t_]=y[t]/.hog[[1]];
r[t_]=r[t]/.hog[[1]];
(Null^2)
ParametricPlot3D[{t, r[t], y[t]}, {t, 0, 20}, Boxed->False, BoxRatios->{5, 2, 2}, AxesLabel->{t, r, y}]

⁃Graphics3D⁃

Do More Practice:  
1- Change the speed of the adjustments in the goods and money markets from .4 and .6 to .2 and .9, respectively.   Re-run th eprogram and observe changes in the solution and the adjustment path.
2- Change the marginal propensity to consume in the consumption function from .9 to .7.  Change the coefficient of income in the demand for money equation from .25 to .30.  Re-run the program and observe changes in the solution and the adjustment path.

Example:  Application to IS-LM Model
Assume that a macro model of an economy is represented by the following equations: c = 350 + .8y,  i = 75 - 3r, y = c + i,  Md = .25y - .8r, and Md = Ms = 490, where the variables and the relationships are defined as in the previous example.    Suppose that the speed of the adjustments in the goods market and money market are proportional to the excess demand in the corresponding markets, with a proportionality factor of .4 for the goods market and a proportionality factor of .6 for the money markets.  Given that the initial values of y and r are y(0) = 1500 and r(0) = 9, find the time path of y and r and analyze the stability of the model.  
Solution:  The equilibrium condition in the goods market requires y = c + i = 425 + .8y - 3r and the equilibrium condition in the money market requires .25y - .8 r = 490.  The slope of the IS curve is, (dr/dy)=  .(2/(-3)) = - .066 and the slope of the LM curve is (dr/dy)= \(TraditionalForm\` .25\/ .8\)= .312.   IS and LM have the conventional slopes and the trace condition for the stability of the model is satisfied.  The determinant condition for the stability of the model requires that the slope of the LM curve exceeds the slope of the IS curve.  This condition is also satisfied.   Therefore, the model is a stable model.  However, the solution of the system of differential equations using Mathematica shows that the solution involves imaginary numbers.  Therefore, the time path of variables spirals to equilibrium.  To solve the problem using Mathematica and to observe the time path to equilibrium, use the following commands,
ClearAll[y,c,i, r,Md, Ms, eel]
c=350+.8y[t]; i=75-3r[t];  Md = .25y[t]-.8r[t]; Ms = 490;
eel=DSolve[{y'[t]==.4(c+i-y[t]), r'[t]==.6(Md -Ms), r[0]==9, y[0]==1500}, {y[t], r[t]}, t]
\({{y[
              t] -> \((\(\(0.`\)\(\)\) - 
                  0.021978021978021976`\ i)\)\ \((\((\(\(5941.751930197019`\)\(\)\) - 
                        11125.`\ i)\)\ e\^\(\((\(-0.28`\) - 0.37416573867739417`\ i)\)\ t\) - \((\(\(5941.751930197019`\)\(\)\) + 
                        11125.`\ i)\)\ e\^\(\((\(-0.28`\) + 0.37416573867739417`\ i)\)\ t\) + \((\(\(0.`\)\(\)\) + 
                        90500.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\), 
          r[t] -> \((\(\(0.`\)\(\)\) - 
                  0.016483516483516484`\ i)\)\ \((\((\(\(5945.493587583793`\)\(\)\) - 
                        2.`\ i)\)\ e\^\(\((\(-0.28`\) - 0.37416573867739417`\ i)\)\ t\) - \((\(\(5945.493587583793`\)\(\)\) + 
                        2.`\ i)\)\ e\^\(\((\(-0.28`\) + 0.37416573867739417`\ i)\)\ t\) + \((\(\(0.`\)\(\)\) + 
                        550.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\)}}\)
the following commands will yield the time path of y(t)
Short[y[t_]=y[t]/.eel[[1]]]
Plot[y[t], {t, 0, 5}]
(0.-0.021978 i) («1»)

⁃Graphics⁃
the following commands will yield the time path of r(t)
Short[r[t_]=r[t]/.eel[[1]]]
Plot[r[t], {t, 0, 10}]
\(\((\(\(0.`\)\(\)\) - 
              0.016483516483516484`\ i)\)\ \((\((\(\(5945.493587583793`\)\(\)\) - 
                    2.`\ i)\)\ e\^«1» - «1» + \((\(\(0.`\)\(\)\) + 
                    550.`\ i)\)\ e\^«1»)\)\)

⁃Graphics⁃
ParametricPlot3D[{t, r[t], y[t]}, {t, 0, 10}, Boxed->False, BoxRatios->{5, 2, 2}, AxesLabel->{t, r, y}]

⁃Graphics3D⁃
To graph the IS-LM on the standard two dimensional space, use the following commands:
To solve for the IS equation from the equilibrium condition in the goods market,
ise=Solve[y==425+.8y-3r, r]
{{r->-0.333333 (-425.+0.2 y)}}
To plot the IS,
is[y_]=r/.ise[[1]];
g1=Plot[is[y],{y, 1800, 2200}]	

⁃Graphics⁃
To solve for the LM equation from the equilibrium condition in the money market,
lme=Solve[490==.25y-.8r, r]
{{r->-1.25 (490.-0.25 y)}}
To graph the LM,
lm[y_]=r/.lme[[1]]
g2 =Plot[lm[y], {y, 1800, 2200}]
-1.25 (490.-0.25 y)

⁃Graphics⁃
To graph the parametric plot of IS and Lm over t on the two dimensional y - r space,
g3=ParametricPlot[{y[t], r[t]}, {t, 0, 20}]


⁃Graphics⁃
To show the three graphs together on the y - r space, 
Show[g1, g2, g3]

⁃Graphics⁃

Do More Practice:  
1- Change the speed of the adjustments in the goods and money markets from .4 and .6 to .2 and .9, respectively.   Re-run the program and observe changes in the solution and the adjustment path.
2- Change the marginal propensity to consume in the consumption function from .8 to .7.  Change the coefficient of income in the demand for money equation from .25 to .30.  Re-run the program and observe changes in the solution and the adjustment path.

Example:
Solve the following system of linear homogeneous differential equations and analyze the stability of the system:  y'1 = y1 + .8y3;  y'2 = y1 -.2y2;  y'3 = -.3y3.
Solution:  The system can be written as Y' = AY = (({{-1, 0, .8}, {1, - .2, 0}, {0, 0, - .3}}))Y.  The eigenvalues of the A matrix of the system are,
Eigenvalues[(({{-1, 0, 0}, {1, - .2, 0}, {0, 0, - .3}}))]
{-1.,-0.3,-0.2}
the eigenvalues are all real and negative.  Therefore, the system will be a stable one.  To solve the system and to graph the time path of the variables,
 Clear[y1,y2,y3,ray]
ray=DSolve[{y1'[t]==-y1[t]+.8y3[t], y2'[t]==y1[t]-.2y2[t], y3'[t]==-.3y3[t], 
      y1[0]==y2[0]==y3[0]==1}, {y1[t], y2[t], y3[t]}, t]
\({{Null\ \((y1[t] -> 
                0.14285714285714285`\ e\^\(\(-1.`\)\ t\)\ \((\(-1.`\) + 8.`\ e\^\(0.7`\ t\))\))\), 
          Null\ \((y2[t] -> 
                0.03571428571428571`\ e\^\(\(-1.`\)\ t\)\ \((\(\(5.`\)\(\)\) - 320.`\ e\^\(0.7`\ t\) + 
                      343.`\ e\^\(0.8`\ t\))\))\), 
          Null\ \((y3[t] -> e\^\(\(-0.3`\)\ t\))\)}}\)
y1[t_]=y1[t]/.ray[[1]]
	Plot[y1[t], {t, 0, 10}]
0.14285714285714285` (e^(-1.` t)) (-1.` + 
            8.` (e^(0.7` t)))

⁃Graphics⁃
y2[t_]=y2[t]/.ray[[1]]
Plot[y2[t], {t, 0, 10}]
\(0.03571428571428571`\ e\^\(\(-1.`\)\ t\)\ \((\(\(5.`\)\(\)\) - 320.`\ e\^\(0.7`\ t\) + 
            343.`\ e\^\(0.8`\ t\))\)\)

⁃Graphics⁃
y3[t_]=y3[t]/.ray[[1]]
Plot[y3[t], {t, 0, 10}]
(e^(-0.3` t))

⁃Graphics⁃
ParametricPlot3D[{y1[t], y2[t], y3[t]}, {t, 0, 20}, Boxed->False, BoxRatios->{5, 2, 2}]

⁃Graphics3D⁃

Do More Practice:  
1- Change the A matrix of the system to (({{1, 0, .8}, {1, .2, 0}, {0, 0, .3}}))and re-run the program.   Compare the results to those in the text.

2-  Change the A matrix of the system to (({{1, 0, .8}, {1, - .2, 0}, {0, 0, .3}}))and re-run the program.   Compare the results to those in the text.

3-  Change the A matrix of the system to (({{1, 0, .8}, {1, .2, 0}, {1, 0, .3}}))and re-run the program.   Compare the results to those in the text.
4- Compare the effects of the changes in the eigenvalues of a system on the stability of the system and the time path of the variables using the practices above.

 Problems:
1.  Analyze the stability of the following differential equations
	a)	y' + 2y = 0
	b)	y' = 2y + 4
	c)	y" + 2y' + y = e
	d)	y" - 2y' + y = (e^t)			       
	e)	y" ' + 6y" + 11y' + 6y = (e^t)
	f)	y" ' + 2y" - y' - 2y = 0
	g)	y" ' - y' + y = 0	
2.  Analyze the stability of the markets for the following models
	a)	Qd = 100 -2p
		Qs = -20 + 3p		
	b)	Qd = 100 - 2p
		Qs = -20 + p
	c)	Qd = 100 - 4p
		Qs =   40  -  p			
	d) 	Qd = 100 -   p
		Qs = 160 - 4p	
3.  A simple macro model of an economy is represented by the following equations: C = 150 + .9y,  I = 100 - 4r, y = C + I,  Md = .2y  - 1.2r, and Md = Ms = 400.  The first equation expresses consumption (C) as a function of income (y) ; the second equation expresses investment ( I ) as a function of interest rate (r); the third equation is the equilibrium condition in the goods market; the fourth equation expresses demand for money (Md) as a function of income and interest rate; and the last equation is the equilibrium condition in the money market.  Suppose that the speed of the adjustments in the goods market and money market are proportional to the excess demand in the corresponding markets, with a proportionality factor of .3 for the goods market and a proportionality factor of .8 for the money markets.  Given that the initial values of y and r are y(0) = 1500 and r(0) = 10, find the time path of y and r and analyze the stability of the model.  
4- Solve and analyze the stability of the following system of linear homogeneous differential equations:  y1' = y2 - .8y3;  y2' = y1 + y3;  y3' = y2 + y3.
5-  Solve and analyze the stability of the following system of linear homogeneous differential equations:  y1' =2y1 - y2 + .8y3;  y2' = -.6y1 + y3;  y3' = 1.1y2 - .3y3.	
	
Appendix 1:  Existence of Definite Solution
One very important theory of the system of differential equations is concerned with the variety of solutions and the general solution to the system of differential equations.  According to this theory, if y1(t) and y2(t) are soultions of y'(t) = Ay(t), then for any arbitrary constants c1 and c2,  y(t) = c1y1(t)  + c2y2(t) is also a solution.  This equation is called the general solution to the system of differential equations.

Example:  Solve the system of equations in the preceding example for the initial values, y(o) = (({{1}, {9}})).
Solution:  The general solution to the probelm is,  (({{y1 (t)}, {y2 (t)}})) =c1 (e^(-t))(({{-1}, {3}})) +  c2 (e^(3  t))(({{1}, {1}})).  To find the solution of the initial-value problem, set t = 0 in the general solution, set (({{y1 (0)}, {y2 (0)}}))= (({{1}, {9}})), and solve for c1 and c2.  That is, (({{1}, {9}}))= (({{-1, 1}, {3, 1}})).(({{c1}, {c2}})).  The solution to this system of equations is,
Solve[(({{-1, 1}, {3, 1}})).(({{c1}, {c2}}))==(({{1}, {9}})),  {c1,  c2}]
{{c1->2,c2->3}}
This results in the initial-value solutions of y1(t) = -2(e^(-t)) + 2(e^(3  t)) and y2(t) = 6(e^(-t)) + 3(e^(3  t)).  
To solve the same problem using Mathematica, type
ClearAll[y1, y2, cat]
cat=DSolve[{y1'[t]==2y1[t]+y2[t], y2'[t]==3y1[t], y1[0]==1, y2[0]==9}, {y1[t], y2[t]}, t]
\({{y1[
              t] -> e\^\(-t\)\ \((\(-2\) + 
                  3\ e\^\(4\ t\))\), 
          y2[t] -> 
            3\ e\^\(-t\)\ \((2 + e\^\(4\ t\))\)}}\)
Once the problem is solved, use the replacement rule "/." to define the function y1[t] and y2[t] and graph the time path of the variables.
y1[t_]=y1[t]/.cat[[1]]
y2[t_]=y2[t]/.cat[[1]]
ParametricPlot3D[{y1[t], y2[t], t}, {t, 0, 10}, Boxed->False, BoxRatios->{4,2,2}, AxesLabel->{y1,y2,t}]
(e^(-t)) (-2 + 
            3 (e^(4 t)))
3 (e^(-t)) (2 + (e^(4 t)))

⁃Graphics3D⁃

Do More Practice:  
1- Change the coefficients of the A matrix in the preceding example to (({{1, 0}, {1, 4}})), find the eigenvalues of the A matrix, and re-run the program.  Obsrerve the changes in the eigenvalues and their relations to the time path of the solution.
2- Change the coefficients of the A matrix in the preceding example to (({{0, 1}, {3, 4}})), find the eigenvalues of the A matrix, and re-run the program.  Obsrerve the changes in the eigenvalues and their relations to the time path of the solution.
3- Change the coefficients of the A matrix in the preceding example to (({{1, -1}, {4, 4}})), find the eigenvalues of the A matrix, and re-run the program.  Observe the changes in the eigenvalues and their relations to the time path of the solution.and redo the example.

A system of n differential equations has a distinct solution for the given initial values if the matrix A of the system has n distinct eigenvalues.  If the system does not have enough distinct eigenvalues and eignevectors, the variety of solutions is reduced and the initial-value problem may not be solvable.  However, we can search for other solutions that will provide us with a definite solution to the initial value problem.  This point will be illustrated through the following example.

Example:  Solve the system of differential equations  y'(t) = (({{1, 3}, {0, 1}})) y (t) for the initial values y(0) = (({{1}, {1}})).
Solution:  The eigenvalues and eigenvectors of the A matrix are, 
Eigensystem[(({{1, 3}, {0, 1}}))]
{{1,1},{{1,0},{0,0}}}
Since the A matrix is of order 2x2, to solve the initial-value problem we need two distinct eigenvectors.  However, the matrix A does not have two distinct eigenvalues, which implies it does not have two distinct eigenvectors.   This does not necessarily imply that we can not find two distinct eigenvectors and solve the initial-value proble.  If the matrix A does not have enough distinct eigenvalues, and there is a deficiency in the number of distinct or linearly independent eigenvectors, we have to search for other solutions that will provide us with enough independent solutions to y' (t) = Ay (t) so that it could be solved for the initial values.  One possibility is to look for functions of the form, y (t) = (u t  + v )(e^kt) or  y (t) = (u (t^2) + v t + w)(e^kt),  etc.	  
Let us explore a solution of the form y (t) = (u t  + v )(e^kt).   For y (t) = (u t  + v )(e^kt) to be a solution of  y' (t) = Ay (t), the differenttial of y (t) in y (t) = (u t  + v )(e^kt) should be equal to Ay(t).  That is, y' (t) = u(e^kt) + k(ut + v)(e^kt) = Ay(t) = A(ut(e^kt) + v(e^kt)), or u(e^kt) + kut(e^kt) + kv(e^kt) = Aut(e^kt) + Av(e^kt).  By rearranging the terms we get, (Au - ku)t(e^kt) + (Av - kv - u)(e^kt) = 0.   For this equation to be satisfied, the two equations (A - k I)u = 0 and (A - k I)v = u should be true.   The first equation, (A - k I)u = 0, is the solution that we had before.  The second equation, (A - k I)v = u, provides us with a new solution to the system.
In the example above, the only distinct eigenvalue of the system is k = 1, with the corresponding eigenvector (({{1}, {0}})).   To find another solution we use  (A - k I)v = u.  Since we already have k and u, we can find v from the equation  (A - k I)v = u  by solving  ((({{1, 3}, {0, 1}}))- (({{1, 0}, {0, 1}})))(({{v1}, {v2}}))=(({{1}, {0}})),
Solve[(({{0, 3}, {0, 0}})).(({{v1}, {v2}}))==(({{1}, {0}})), {v1,  v2}]

{{v2 -> (1/3)}}
here, v2 = (1/3) and v1 can be any arbitrary number, sy 1.  Therefore, the two solutions are y1 = (e^t)(({{1}, {0}})) and y2 = (te^t)(({{1}, {0}}))+ (e^t)(({{1}, {(1/3)}}))= (e^t)(({{t + 1}, {(1/3)}})).  
The general solution of the problem is, y(t) = c1y1(t) +c2y2(t) = c1(e^t)(({{1}, {0}})) + c2 (e^t)(({{t + 1}, {(1/3)}})).  Solving this equation for the initial values yields, (({{1}, {1}}))= (({{c1}, {0}}))+ (({{c2}, {(1/3) c2}}))or,  c2 = 3 and c1 = -2.  Therefore. the initial-value problem has a solution y(t) = -2y1(t) + 3y2(t) or (({{y1 (t)}, {y2 (t)}}))= (({{(e^t) + 3  (te^t)}, {(e^t)}})).   Using Mathematica, the solution is,
Clear[y1, y2, ape]
ape=DSolve[{y1'[t] == y1[t]+3y2[t], y2'[t]==y2[t], y1[0]==1, y2[0]==1}, {y1, y2}, t]
\({{y1 -> \((e\^#1\ \((1 + 3\ #1)\) &)\), 
          y2 -> \((e\^#1 &)\)}}\)
To plot the time path of y1[t] and y2[t], use the replacement rule "/." to define the functions y1[t] and y2[t].
y1[t_]=y1[t]/.ape[[1]]
(e^t) (1 + 3 t)
y2[t_]=y2[t]/.ape[[1]]
(e^t)
ParametricPlot3D[{y1[t], y2[t], t}, {t, 0, 10}, Boxed->False, BoxRatios->{4,2,2}, AxesLabel->{y1,y2,t}]

⁃Graphics3D⁃

Example:  Application to Macroeconomics
The change in national debt of an economy (the budget deficit) is a constant proportion .03 of its GNP and the GNP growth of the economy is 8%.  Given the initial GNP is 6000 and the debt is 3000, find the time path of GNP and the debt.
Solution:  Lets call the GNP,  y and the national debt x.  The problem is to solve the system of equations x' = .03y and y' = .08y, for the initial vali\ues y(0) = 6000 and x(0) = 4000.  The system can be written in matrix form as,  (({{x'}, {y'}})) = (({{0, .03}, {0,  .0  .8}}))(({{x}, {y}})).
ClearAll[ x, y, mot]
mot=DSolve[{x'[t]==.08y[t], y'[t]==.03y[t], x[0]==4000, y[0]==6000}, {x[t], y[t]}, t]
{{Null (x[t] -> 
                4000.` (-3.` + 
                      4.` (e^(0.03` t)))), 
          Null (y[t] -> 6000.` (e^(0.03` t)))}}
The  plot of national debt over tine  is,
x[t_]=x[t]/.mot[[1]]
4000.` (-3.` + 4.` (e^(0.03` t)))
Plot[x[t], {t, 0, 10}]

⁃Graphics⁃
the plot of GNP  over time is,
y[t_]=y[t]/.mot[[1]]
Plot[y[t], {t, 0, 10}]
6000.` (e^(0.03` t))

⁃Graphics⁃
the parametric plot of debt and GDP over time is, 
ParametricPlot3D[{x[t], y[t], t}, {t, 0, 200}, Boxed->False, BoxRatios->{4,2,2}, AxesLabel->{y1,y2,t}]

⁃Graphics3D⁃
The extent of the affordability of a country's national debt is measured by the ratio of the debt growth to GNP growth.  Since the governments pay interest on debt, a ratio less than one implies that the annual interest payment on debt is less than GNP.  If the ratio is greater than one, the interest payment on debt may turn greater than GNP and the debt may not be affordable.  In this example, the ratio is \(TraditionalForm\` .03\/ .08\) or .375.  That is for every $1 incerase in GNP, the debt is increasing by 37.5 cents.  

Do More Practice:
1- Change the GNP growth to 8% and the debt growth to 3%.  Re-run the program and interpret the solution.
2- Change the GNP growth to 10% and the debt growth to 1%.  Re-run the program and interpret the solution.

Example:  Solve the following system of differential equations, y'(t) = (({{0, -1, 0}, {0, 0, 2}, {0, 0, 0}}))y(t), for the initial values y(0) = (({{1}, {2}, {3}})).
Solution:  Find the eigenvalues and eigenvectors of the A matrix of the system of equations,
Eigensystem[(({{0, -1, 0}, {0, 0, 2}, {0, 0, 0}}))]
{{0,0,0},{{1,0,0},{0,0,0},{0,0,0}}}
The A matrix of the system has only one distinct eigenvalue k = 0, with the corresponding eigenvector u = (({{1}, {1}, {0}})).  Therefore, one solution is, y1(t) = (({{1}, {0}, {0}})).  To find other solutions, we try a solution of the form,  y2(t) = (ut + v)(e^kt).  As proved in the previous example, for  y2(t) = (ut + v)(e^kt) to be a solution the equation (A - k I)v = u should be satisfied.  Since k = 0, (({{0, -1, 0}, {0, 0, 2}, {0, 0, 0}}))(({{v1}, {v2}, {v3}})) = (({{1}, {0}, {0}})).  This equation results in the solution set, v2 = -1, v3 = 0, and v1 any arbitrary number, say v1 = 1.  That is, v = (({{1}, {-1}, {0}})).  Therefore, the second solution of the system of equations is , y2(t) = (ut + v)(e^kt) = (({{1}, {1}, {0}}))t + (({{1}, {-1}, {0}})) = (({{t + 1}, {t - 1}, {0}})).
The system still cannot be solved for the initial values and requires one more solution.   Therefore, we search for a solution of the form,  y'3(t) = (u(t^2) + vt + w)(e^kt).  Differentiating this equation and setting the differential equal to Ay3(t) yields, 	(2ut + v)(e^kt) + k(u(t^2) + vt + w)(e^kt) =  A(u(t^2) + vt + w)(e^kt).  By rearranging the terms we get, (Au - ku)(t^2)(e^kt) + (Av - kv - 2u)t(e^kt) + (Aw - kw - v)(e^kt) = 0.   For this equation to be zero the following three equations should satisfy:
	Au - ku = 0,   or   (A - kI)u = 0
	Av - kv = 2u,  or  (A - kI)v = 2u
	Aw -kw = v,   or  (A - kI)w = v
The first equation results in the eigenvalue k = 0, with the corresponding eigenvector u = (({{1}, {0}, {0}})).  The second equation yields, (({{0, -1, 0}, {0, 0, 2}, {0, 0, 0}}))(({{v1}, {v2}, {v3}}))= (({{2}, {0}, {0}})), or  v = (({{v1}, {v2}, {v3}}))= (({{1}, {-2}, {0}})).  The third equation yields (({{0, -1, 0}, {0, 0, 2}, {0, 0, 0}}))(({{w1}, {w2}, {w3}}))= (({{1}, {-2}, {0}})), or w = (({{w1}, {w2}, {w3}}))= (({{1}, {-1}, {-1}})).  Therefore, the third solution is y3(t) = (({{1}, {0}, {0}}))(t^2)+ (({{1}, {-2}, {0}}))t + (({{1}, {-1}, {-1}})).  
 Now that we have found three linearly independent solutions to the system, we can write the general solution as, 	y(t) = c1y1(t) + c2y2(t) + c3y3(t) = c1(({{1}, {0}, {0}}))+ c2((({{1}, {0}, {0}}))t +  (({{1}, {-1}, {0}})))+ c3((({{1}, {0}, {0}}))(t^2) + (({{1}, {-2}, {0}}))t + (({{1}, {-1}, {-1}}))).  Solving this equation for the initial values y(0) = (({{y1 (0)}, {y2 (0)}, {y3 (0)}})) = (({{1}, {2}, {3}})) yields c1 = 3,  c2 = 1, and c3 = -3.   Therefore, the inital-value solution is, y(t) = (({{-3 (t^2) - 2  t + 1}, {6  t  + 2}, {3}})).  Using Mathematica, the solution is,
ClearAll[y1,y2,y3]
DSolve[{y1'[t]  == 0, y2'[t]==-y2[t], y3'[t]==2y3[t], y1[0]==1, y2[0]==2, y3[0]==3}, {y1,y2,y3}, t] 
{{y1 -> (1 &), 
          y2 -> (2 (e^(-#1)) &), 
          y3 -> (3 (e^(2 #1)) &)}}

Appendix 2:  Solution of the Systems of First Order, Nonhomogeneous Linear Differential Equations
The resemblance of the system of first order, nonhomogenous linear differential equations y'(t) = Ay(t) + f(t) to the first order, nonhomogenous linear differential equation y'(t) = ay(t) + f(t) suggests that a similiar solution may be applicable here.  However, before finding solution, two concepts which are very useful in solving systems of nonhomogeneous differential equations will be discussed below. 

Definition 1. Fundamental Matrix:    A nxn matrix, M(t), whose columns are linearly independent solutions of the system y'(t) = Ay(t) at to is called a fundamental matrix of y'(t) = Ay(t) at to.

Example:  Find a fundamental matrix for the following system of differential equations y'(t) = (({{1, -1}, {0, 2}}))y(t).          
Solution:  The system has two distinct eigenvalues k1 = 1 and K2 = 2, with the corresponding eigenvectors u1 = [(({{1}, {0}}))and u2 = (({{1}, {-1}})).  Therefore, the soulutions of the system of equations are, y1(t) = (({{1}, {0}}))(e^t) and  y2(t) = (({{1}, {-1}}))(e^(2  t)).   The solution vectors are linearly independent and the fundamental matrix of the system is  y(t) = (({{(e^t), (e^(2  t))}, {0, -(e^(2  t))}})), (note that,  if the eigenvalues are distinct, the corresponding eigenvectors will be linearly independent.).   Clearly,  the fundamental matrices are not unique.  Any combination of solution vectors which are linearly independent will give a fundamental matrix for the system.  It is also easy to prove that cM(t), where c is a constant vector, is a solution to the system y'(t) =Ay(t).

Definition 2. Normalized Fundamental Matrix:   A matrix formed as Mn(t) = M(t)(M^(-1))(to) is called a normalized fundamental matrix.  From the definition of normalized fundamental matrices, it is easy to see that they are identity matrices at t = 0.

Example:  Find the normalized fundamental matrix for the system of differntial equations of the preceding example.
Solution:  The fundamental matrix of the preceding example is,  	M(t) = (({{(e^t), (e^(2  t))}, {0, -(e^(2  t))}})).  At t = 0, M(0) = (({{1, 1}, {0, -1}})) and,  (M^(-1))(0)  =  (({{1, 1}, {0, -1}})).  The normalized matrix is Mn(t) = M(t)(M^(-1))(0)= (({{(e^t), (e^t) - (e^(2  t))}, {0, (e^(2  t))}})).  This  matrix will be an identity matrix at t = 0.

Definitions 1 and 2 help us to find a solution to the nonhomogenous differential equations.  Associated with every nonhomogenous system is a complementary homogenous system, obtained by setting f(t) = 0.  If M(t) is the fundamental matrix of the complementary system, we search for a vector u(t) such that,  y(t) = M(t)u(t)  is a solution to the nonhomogenous system.  If  y(t) = M(t)u(t) is a solution to the system, its differential will correspond to equation y'(t) = Ay(t) + f(t).  Differential of  y(t) = M(t)u(t) yields, y'(t) = M'(t)u(t) + u'(t)M(t).  Since M(t) is a fundamental matrix of the homogenous system, M'(t) = AM(t) and y'(t) = M'(t)u(t) + u'(t)M(t) can be written as, y'(t) = A(M(t)u(t)) + u'(t)M(t) = Ay(t) + f(t).  This equation implies that, M(t)u'(t) = f(t), or u'(t) = (M^(-1))(t)f(t).  The preceding equation together with the equation  y(t) = M(t)u(t) provide a solution to the system.  However, this is only one of the many particular solutions of the system.   A general solution will be, y(t) = yp(t) + yc(t)C = M(t)u(t) + M(t)C, where C is any vector of constants.   The equation y(t) = M(t)u(t) + M(t)C is the solution of  y'(t) = Ay(t) + f(t) for any choice of C.

Example:  Find the solution of the follwing system of nonhomogeneous differential equations y'(t) = (({{1, -1}, {0, 2}}))y(t) + (({{(e^(2  t))}, {1}})).  
Solution:  To find the fundamental matrix of the complementary homogenous system, as defined in the definition 1 above, find the eigensytem of the A matrix, 
Eigensystem[(({{1, -1}, {0, 2}}))]
{{1,2},{{1,0},{-1,1}}}
the fundamental matrix of the system is  M(t) = (({{(e^t), -(e^(2  t))}, {0,   (e^(2  t))}})).  Use the expression  u'(t) = (M^(-1))(t)f(t) to find u(t) of the solution y(t) = yp(t) + yc(t)C = M(t)u(t) + M(t)C.  That is, u'(t) = (({{(e^(-t)), (e^(-t))}, {0, (e^(-2 t))}}))(({{(e^(2  t))}, {1}}))=(({{(e^t) + (e^(-t))}, {(e^(-2 t))}})).  To find u(t), integrate u'(t) in the interval zero to t, u(t) = (∫_0^t)(({{(e^t) + (e^(-t))}, {(e^(-2 t))}})) dt , 
Integrate[{E^t+E^(-t), E^(-2t)},  {t, 0, t}]
{(e^(-t)) (-1 + (e^(2 t))), 
        (1/2) - ((e^(-2 t))/2)}
The particular solutin of the system is, yp(t) = M(t)u(t) =  (1/2)(({{(e^t), -(e^(2  t))}, {0,   (e^(2  t))}}))(({{2  (e^t) - 2  (e^(-t))}, {1 - (e^(-2 t))}}))= (1/2) ((e^(2  t)) - 1)(({{1}, {1}})).  This solution to the problem is only one of the many particular solutions.  A general solution of the problem is, y(t) = yp(t) + yc(t)C = M(t)u(t) + M(t)C = (1/2)((e^(2  t)) - 1)(({{1}, {1}}))+(({{(e^t), -(e^(2  t))}, {0,    (e^(2  t))}}))(({{c1}, {c2}})) where c1 and c2 are any arbitrary constants.  
To solve the same problem using Mathematica, use the command "DSolve",
ClearAll[y1,y2]
DSolve[{y1'[t]==y1[t]-y2[t] + E^(2t), y2'[t]==2y2[t]+1}, {y1[t], y2[t]}, t]
{{y1[t] -> 
            (1/2) (-1 + 2 (e^(2 t)) + 
                  2 (e^t) C[1] + 2 (e^t) C[2] - 
                  2 (e^(2 t)) C[2]), 
          y2[t] -> 
            (1/2) (-1 + 2 (e^(2 t)) C[2])}}

Example:  Find the solution of the following system of differential equations, y(t) = (({{1, 0, 1}, {0, 0, 1}, {0, 0, 0}}))+(({{2  (e^t)}, {t}, {1}})).
Solution:  The eigensystem of the  A matrix is, 
Eigensystem[(({{1, 0, 1}, {0, 0, 1}, {0, 0, 0}}))]
{{0,0,1},{{0,1,0},{0,0,0},{1,0,0}}}
the matrix has two distinct eigenvalues k1 = 0 and K2 = 1,  with the corresponding eigenvectors u1 = (({{0}, {1}, {0}}))and u2 =(({{1}, {0}, {0}})).  Therefore, the two solutions corresponding to k1 and k2 are, y1(t) = (({{0}, {1}, {0}})) and y2(t) = (({{1}, {0}, {0}}))(e^t).  To find the general solution of the system, we need one more distinct solution.  Let us try a solution of the from y3(t) = (ut + v)(e^t) (see the section on the existence of definite solution).   (A - kI)v = u yields the v vector, v =(({{-1}, {1}, {1}}))and the third solution is, y3(t) = (({{0}, {1}, {0}}))t+ (({{-1}, {1}, {1}}))= (({{-1}, {t + 1}, {1}})).  The fundamental matrix of the system is M(t) = (({{0, (e^t), -1}, {1, 0, t + 1}, {0, 0, 1}}))and (M^(-1))(t) = (({{0, 1, -t - 1}, {(e^(-t)), 0, (e^(-t))}, {0, 0, 1}}))
									
and u'(t) = (M^(-1))(t)f(t) =  (({{0, 1, -t - 1}, {(e^(-t)), 0, (e^(-t))}, {0, 0, 1}})) (({{2  (e^t)}, {t}, {1}}))= (({{-1}, {2 + (e^(-t))}, {1}})).  Integrate u' (t) in the interval zero to t to get u(t), u(t) = (({{-t}, {2  t - (e^(-t))}, {t}})+1).  Using this equation we can find y(t) as,
y(t) = M(t)u(t) = (({{0, (e^t), -1}, {1, 0, t + 1}, {0, 0, 1}}))(({{-t}, {2  t - (e^(-t))}, {t}})+1)= (({{2  (te^t) + (e^t) - t - 1}, {(t^2)}, {t}})).  This, is only one of the many particular solutions of the system.  To find the general solution, y(t) = yp(t) + yc(t)C = M(t)u(t) + M(t)C, where C is any arbitrary constant.  If the initial value of y(to) = yo is also given, the initial value problem y'(t) = Ay(t),  y(to) = yo has the unique solution y(t) = yp(t) + Mn(t)yo,  where Mn(t) is the normalized fundamental matrix of the complementary system.

Example:  Solve the system of differential equations y'(t) = (({{1, -1}, {0, 2}}))y(t) + (({{(e^(2  t))}, {1}})) for the initial value, y(0) = (({{1}, {1}})).
Solution:  The solution of the complementary system of the problem is y(t) = (1/2)((e^(2  t)) -1)(({{1}, {1}})).  Normalize the fundamental matrix of the system to get Mn(t),  Mn(t) = M(t)(M^(-1))(0)= (({{(e^t), (e^(2  t))}, {0, -(e^(2  t))}}))(({{1, 1}, {0, -1}}))= (({{(e^t), (e^t) - (e^(2  t))}, {0, (e^(2  t))}})).  Therefore, the solution to the initial-value problem is y(t) = yp(t) + Mn(t)yo = (1/2) ((e^(2  t))  - 1)(({{1}, {1}}))+ (({{(e^t), (e^t) - (e^(2  t))}, {0, (e^(2  t))}}))(({{1}, {1}}))= (1/2)(({{4  (e^t) - (e^(2  t)) - 1}, {3  (e^(2  t)) - 1}})).
To solve the same problem using Mathematica, 
ClearAll[y1,y2]
DSolve[{y1'[t]==y1[t]-y2[t] +E^(2t), y2'[t]==2y2[t] +1, y1[0]==1, y2[0]==1}, {y1[t], y2[t]}, t]
{{y1[t] -> 
            (1/2) (-1 + 
                  4 (e^t) - (e^(2 t))), 
          y2[t] -> 
            (1/2) (-1 + 3 (e^(2 t)))}}

Example:  Application to Markets
One of the major characteristics of a perfectly competitive industry is that the firms have easy entry or exit to the industry.  Suppose that the firms' entry or exit in a perfectly competitive industry (the change in the number of firms, (dn/dt)) is proportional to the gap between the market price (p) and the minimum average total cost (atc), (dn/dt) = n' = λ(p - atc), λ > 0.  That is, if p - atc > 0, the positive economic profit attracts more entry to the industry and results in lowers.  The opposite occurs when p - atc < 0.  In addition, assume that changes in the market price is proportional to excess demand in market, (dP/dt) = p' = μ(Qd - Qs),  μ > 0, where Qd and Qs represent market demand and market supply, respectively.  λ and μ represent the speed of adjustment for the number of firms and price, respectively.  Given that Qd = 100 - 2p, Qs =  5n, atc = 20, λ = .1 and μ = .2 formulate the price behaviour in the market, solve the resulting system of differential equations, and plot the time path of the price.  The initial values of p and n are given as po = 25 and no = 15.  Note that, in this problem, it is assumed that the firms in the industry are of the uniform size and that Qs depends on the number of firms in the industry.
Solution:  The system of the differential equations can be expressed as p' = .2(100 -2p -5n) and n' = .1(p - 20).  In matrix form the system can be written as (({{p'}, {n'}})) = (({{- .4, -1}, {.1, 0}}))(({{p}, {n}})) + (({{20}, {-2}})).  The steady-state value of p and n are obtained by setting p' = n' = 0.  Using Mathematica,
ClearAll[Qd,  Qs,  p,  n,  chimp ]
Solve[(({{- .4, -1}, {.1, 0}})).(({{p}, {n}}))+(({{20}, {-2}}))==(({{0}, {0}})), {p,  n}]
{{Null (p->20.),Null (n->12.)}}
the steady-state values of p and n are 20 and 12, respectively.  To solve the system of differential equations using Mathematica,
Qd = 100 - 2p[t]; Qs = 5n[t];
chimp = DSolve[{p'[t]==.2(Qd-Qs), n'[t]==10(p[t]-20), p[0]==25, n[0]==15}, {p[t], n[t]}, t]
\({{Null\ \((p[
                  t] -> \((\(\(0.`\)\(\)\) - 
                      0.010040160642570285`\ i)\)\ \((\((\(\(63.118935352238`\)\(\)\) + 
                            249.`\ i)\)\ e\^\(\((\(-0.2`\) - 3.1559467676119`\ i)\)\ t\) - \((\(\(63.118935352238`\)\(\)\) - 
                            249.`\ i)\)\ e\^\(\((\(-0.2`\) + 3.1559467676119`\ i)\)\ t\) + \((\(\(0.`\)\(\)\) + 
                            1992.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\))\), 
          Null\ \((n[
                  t] -> \((\(\(0.`\)\(\)\) + 
                      0.002008032128514057`\ i)\)\ \((\((\(\(3992.2726610290533`\)\(\)\) - 
                            747.`\ i)\)\ e\^\(\((\(-0.2`\) - 3.1559467676119`\ i)\)\ t\) - \((\(\(3992.2726610290533`\)\(\)\) + 
                            747.`\ i)\)\ e\^\(\((\(-0.2`\) + 3.1559467676119`\ i)\)\ t\) - \((\(\(0.`\)\(\)\) + 
                            5976.`\ i)\)\ e\^\(\((\(\(0.`\)\(\)\) + 0.`\ i)\)\ t\))\))\)}}\)
p[t_]=p[t]/.chimp[[1]];
Plot[p[t], {t, 0, 10}]

⁃Graphics⁃