<< Graphics`

Using Mathematica to Solve Differential Equations

Single differential equations

Symbolic solutions

The basic command in Mathematica for solving a differential equation symbolically is DSolve. This works for single differential equations to find a general solution, initial value problems, and systems of differential equations. The one thing you must be careful to do when working with differential equations in Mathematica is to explicitly specify each function and its independent variable.  So, in order to enter the differential equation:  y'' - 2 y' + y = sin x, you would have to enter it as: y''[x] - 2 y'[x] + y[x] == 
        Sin[x] (notice that I also converted it into "Mathematica-speak" at the same time; you must you == instead of =, square brackets, and capitalize all built-in functions).

General solutions

DSolve[y''[x] + 2 y'[x] + y[x] == Sin[x], y[x], 
      x]

Notice that, in addition to specifying the differential equation to solve, you must also specify the function to solve for (y[x]) and the independent variable (x).  In the solution, C[1] and C[2] are the two arbitrary constants for the general solution to this differential equation.

If all you need to do is see your answer, the above works fine.  However, normally, you want to take this solution and do other things with it (plug numbers into it, graph it, make it into small table decorations, etc.).  The form that Mathematica gives you the answer in here (y[x] ->   ... ) is what is called a "replacement rule".  On the one hand, this is a very powerful form; you can use it in some very fancy ways.  However, it is a little tricky to use if you aren't pretty familiar with Mathematica.  So, most of the time we will want to simply assign this solution to a function we can work with later:

Clear[ySol, diffEqn]

(You should get in the habit of Clear-ing all new names that you are about to use.  Otherwise, old definitions you have forgotten about could really screw things up...)

diffEqn = y''[x] + 2 y'[x] + y[x] == Sin[x]

ySol[x_] = 
      y[x] /. First[
          DSolve[y''[x] + 2 y'[x] + y[x] == Sin[x], y[x], 
            x]]

This defines a function ySol[x] to be our solution function. (The /. is an operator that says "replace all the y[x]'s to the left of this using the rule that follows", which in this case is the solution to the differential equation, which said y[
          x] -> (e^(-x)) C[
              1] + (e^(-x)) x C[2] - Cos[x]/2.  The First command strips out the extra {}'s in the solution.) 

Now, we check our answer:

diffEqn /. y -> ySol // FullSimplify

The above line says "Take the original diffEqn, replace all the y's in it with ySol (our solution) and then simplify the result." If the answer is "True" then the solution satisfies the differential equation.

*Useful hint: If you are getting complex solutions (which you presumably don't want), read this section.

Try adding one of the following to the end of the DSolve line:
//ComplexExpand//Simplify
OR
//ComplexExpand//FullSimplify

Unfortunately, even this doesn't always work.At that point,you have 4 options (that I can think of, anyway):

1) Solve the problem numerically using NDSolve (probably the best idea, if you can get enough accuracy).
2) Plot the graph anyway (it will just have some blank spots in it)
3) Assume the imaginary part is ignorable and take the Real part of your answer (using the Re function). To be honest,since I'm not entirely sure howMathematica's algorithms work for diff equations, I would be a little hesitant to do this without knowing a bit more. (Notice that this is really a fancier version of option 2.)
4) Use some other method to solve the problem symbolically (i.e., solve the characteristic polynomial numerically, then use what you know about the form of the homogeneous solution and the particular solution). This is probably the best way, but it can be a lot of extra work (and probably won't get you much extra payoff). I would save this for when you are going to do some serious work with that particular diff. equation.

Initial value problems

Most of the time when you work with a differential equation on a computer system, you will be interested in "initial value problems", rather than just finding the general solution. This is handled very simply: you just include the initial conditions in a list with your differential equation to solve. So, to solve y'' - 2 y' + y = sin x, y(0) = 1,  y' (0) = -2, you would:

Clear[x, y, ySol, diffEqn]

diffEqn = y''[x] + 2 y'[x] + y[x] == Sin[x]

ySol[x_] = 
      y[x] /. First[
          DSolve[{y''[x] + 2 y'[x] + y[x] == Sin[x], 
              y[0] == 1, y'[0] == -2}, y[x], 
            x]]

(Again, be sure to use == instead of = everywhere inside DSolve.)

We could graph this solution curve:

Plot[ySol[x], {x, 0, 4  π}]

What if we wanted to graph several solution curves, for different initial conditions?  Let's examine the impact of curves with different initial y values.  To allow ourselves the most flexibility, let's use some dummy variables:

Clear[y0, v0]

ySol[x_] = 
      y[x] /. First[
          DSolve[{y''[x] + 2 y'[x] + y[x] == Sin[x], 
              y[0] == y0, y'[0] == v0}, y[x], x]]

Plot[
      ySol[x] /. {y0 -> -1, v0 -> 2}, {x, 0, 
        4  π}]

This says to plot ySol again, but replace all the y0's with -1 and v0's with 2.  A more interesting question would be to hold the v0's constant (say at v0 = 1) and let the y0's vary from -5 to 5.  There is a really cool way to do this:

Clear[v0, y0]

v0 = 1; y0 = Range[-5, 5, 1]

Plot[Evaluate[ySol[x]], {x, 0, 4  π}, 
      PlotRange -> All]

This plots ySol for all the different values we specified for y0. (You must use the Evaluate command around ySol[x] or Mathematica won't plug in the numbers for y0 before trying to graph and it won't work.)  So, what happens if we let y0 be constant (say y0 = 2) and let v0 vary from -5 to 5?

Clear[v0, y0]

y0 = 2; v0 = Range[-5, 5, 1]

Plot[Evaluate[ySol[x]], {x, 0, 4  π}, 
      PlotRange -> All]

Unknown parameters

This trick is also very handy if you want to examine the effect of different parameters on a differential equation.  So, let's consider a new equation (oh boy...):  x'' + k x = sin t

Clear[diffEqn, x0, v0, k, xSol, x, t]

diffEqn = x''[t] + k x[t] == Sin[t]

xSol[t_] = 
      FullSimplify[
        x[t] /. First[
            DSolve[{diffEqn, x[0] == x0, x'[0] == v0}, x[t], t]], 
        k ∈ Reals]

(The FullSimplify[...,k ∈ 
            Reals] is a useful trick to let Mathematica know that your k is a real number and simplify accordingly.  It doesn't actually matter in this particular problem, but sometimes it does...)

x0 = 1; v0 = 1; k = Range[0.1, 3.1, 0.5]

Plot[Evaluate[xSol[t]], {t, 0, 4  π}]

Notice that I avoid k = 0.  So, what if k is negative?

x0 = 1; v0 = 1; k = Range[-3.1, -0.1, 0.5]

Plot[Evaluate[xSol[t]], {t, 0, 4  π}]

Numeric approxiamtions

Sometimes, Mathematica can't solve a differential equation (if it is evil enough). Or, it might take a very long time for it to solve and you might not really have any need for a complete symbolic solution.  Often, a good numerical approximation is all you really need.  In Mathematica, you you NDSolve (Numeric DSolve) to compute these.  One of the nice things about the way Mathematica does this is that it doesn't just give you a numeric resulet (i.e., x[2] = 3.2222152), it gives you a numerical "interpolating function" which you can then use to compute the value of the function at any point in your computed domain.  Thus, you can graph this solution, find its intercepts, etc.  Not surprisingly, you must provide enough initial conditions to make sure that the solution is unique (i.e., no unknown constants).

Let's solve the following differential equation:  (y^((3))) + y'' + y' = -(y^3) (remember that (y^((3))) means the 3rd derivative), where y(0) = 1,  y' (0) = 0,  y'' (0) = 0.  We must also specify a domain over which we want the solution to be valid, say 0 ≤ x ≤ 20.

Clear[x, y, ySol]

ySol[x_]=y[x]/.First[NDSolve[{Derivative[3][y][x]+y''[x]+y'[x]==-y[x]^3,y[0] == 1,y'[0]==y''[0]==0},y[x],{x, 0, 20}]]

ySol[20]

Plot[ySol[x], {x, 0, 20}]

If you need greater accuracy in your approximation, you can set the WorkingPrecision to be higher (actually, if you really want to investigate this in more detail, you should read up on WorkingPrecision, PrecisionGoal, and AccuracyGoal in the online help; it can get a bit complicated - setting WorkingPrecision to a large number will often help, however).  The following does all internal computations to an internal accuracy of 30 digits (so your answer is probably accurate to about 20 digits):

ySol[x_]=y[x]/.First[NDSolve[{Derivative[3][y][x]+y''[x]+y'[x]==-y[x]^3,y[0] == 1,y'[0]==y''[0]==0},y[x],{x, 0, 20},WorkingPrecision -> 30]]

It's a bit trickier here if you want to graph one of these for a variety of initial conditions:

ySolutionList=Table[y[x]/.First[NDSolve[{Derivative[3][y][x]+y''[x]+y'[x]==-y[x]^3,y[0] == y0,y'[0]==y''[0]==0},y[x],{x, 0, 20}]],{y0, 0, 1.25, 0.25}]

Plot[Evaluate[ySolutionList], {x, 0, 20}]

If you want to solve for multiple initial conditions (like above) or for different values of a parameter, the important thing to remember is that NDSolve will not work with unknown constants in it (unlike DSolve), so you have to plug those constants in multiple times (using the Table command like I did above is a good way to do this).  (You could actually approach DSolve the same way if you wish, but if the differential equation takes a long time to solve, it could really add up to a long solution time in that case).

Plotting solutions and slope fields

Slope fields

To graph a slope-field in Mathematica, you use the command PlotVectorField.  The "vector field" you are actually plotting has the coordinates {1, y'} at each point (since each vector goes "over" 1 and "up" by the slope).  So, to graph the slope field for:  y' = (x/(2 + sin y))

<< Graphics`

PlotVectorField[{1, (x/(2 + Sin[y]))}, {x, -10, 10}, {y, -10, 
        10}]

To clean it up a bit and make it look nicer, try the following options:

PlotVectorField[{1, (x/(2 + Sin[y]))}, {x, -10, 10}, {y, -10, 
        10}, Axes -> True]

If you prefer to make all of the vectors have the same length, try the following:

slopeField = 
      PlotVectorField[{1, (x/(2 + Sin[y]))}, {x, -10, 10}, {y, -10, 
          10}, Axes -> True, ScaleFunction -> (0.4 &), 
        ScaleFactor -> 0.9]

It's a bit complicated how the ScaleFunction -> (0.4 &), 
      ScaleFactor -> 0.9 options work, but if you play around with those numbers a bit you can get something that looks quite nice (be sure that you copy ScaleFunction exactly and then play with the number only; the & is very important).

Putting the solutions and the slope fields together

Often, you will want to plot the solution curves for different initial conditions along with the slope field together on one graph.  You can do that by combining the solution graphing we did above with the slope fields.

Because I picked a pretty nasty equation, DSolve doesn't work like I want, so we will do this one using NDSolve. Let's solve for the different initial values for the curves we want to plot (say from y(0) = -8 to y(0) = 8, every 2 units):

Clear[x, y, y0, ySolutions, plotSolutions]

ySolutions = 
      Table[y[x] /. 
          NDSolve[{y'[x] == (x/(2 + Sin[y[x]])), 
              y[0] == y0}, y[x], {x, -10, 10}], {y0, -8, 8, 
          2}]

plotSolutions = 
      Plot[Evaluate[ySolutions], {x, -10, 10}, 
        PlotStyle -> Blue]

Now, let's combine the slope field with the solution curves:

Show[slopeField, plotSolutions]

This is, of course, hideous.  To fix that, we can explicitly tell it what region to plot:

Show[slopeField, plotSolutions, 
      PlotRange -> {{-10, 10}, {-10, 10}}]

Systems of differential equations

You will notice that I use several different methods to solve systems of equations here.  This is partly because some of them (the sections on eigensystems and exponentiation of operators) are techniques covered in our text; I show you how to work those methods in Mathematica.  Honestly, however, most of the time you will probably use DSolve and NDSolve.  Be sure to read the first section on setting up systems in Mathematica, however.

Setting up the systems - working with matrices

There are a lot of different ways to handle systems of differential equations in Mathematica.  The most simplistic method is to just enter them as lists of differential equations:

DSolve[{x'[t] == x[t] - 2  y[t], y'[t] == 
          2  x[t] + y[t]}, {x[t], y[t]}, t]

This method works (and is easy to understand), but starts to get kind of cumbersome for larger systems (and certainly doesn't win any cool points...).  Personally, I prefer to set systems up from the outset as matrix equations.

To do this, you need to know a few things about how Mathematica handles matrices and vectors.  As far as Mathematica is concerned, a vector is just a list of items: {2, x, y}  A matrix is just a list of lists:  {{1, 2, 3}, {3, 5, 2}, {7, 2, 0}}, where each inside list is a row of the matrix. So this last matrix would look like:

{{1, 2, 3}, {3, 5, 2}, {7, 2, 0}} // MatrixForm

In fact, you can enter matrices directly (so they actually look like matrices) by using the "Create Table/Matrix/Palette" command under the "Input" menu (then just choose "Matrix" and the number of rows/columns) or, if you need a 2 by 2 matrix, there is a button on the standard tool palette for this.

To multiply two matrices, you must use . between them; if you use implied multiplation (or the standard × symbol), Mathematica  won't use matrix multiplication.  So:

(({{2, 4}, {6, 8}})).(({{-3, 4}, {-5, 2}}))//MatrixForm

Or:

A=(({{2, 4}, {6, 8}}));B=(({{-3, 4}, {-5, 2}}))

A . B // MatrixForm

(Sticking //MatrixForm at the end of the line just makes the answer look nice. It isn't necessary.)

One other thing Mathematica does that is kind of nice is that it handles multiplication of a matrix times a vector properly (if you use .), even if the vector isn't really the right kind of vector (i.e., you can be sloppy about row vectors vs. column vectors).  We will take advantage of this to simplify some work that Mathematica would otherwise force you to do.  (I will treat all vectors as row vectors, even when it would be more correct to use column vectors.  Trust me, this simplifies things a lot...)

In order for vector equations to work correctly in Mathematica, you must first execute the following command (this is important, so pay attention):

SetAttributes[Equal, Listable]

(If you are interested, this makes it so the == sign is "listable", i.e., it works with matrices and other lists.)

Now, let's pull all this together and start using matrix notation.  First, let's just set up a system of algebraic equations.  The system:  
({{2  x - 3  y + 4  z = 7}, {2  x - 5  z = 4}, {7  y + 4  x = 9}}) can be written in matrix form with the following definitions:

mA=(({{2, -3, 4}, {2, 0, -5}, {0, 7, 4}}));vC = {7, 4, 9};vX = {x, y, z}

One note about notation:  For clarity sake (and to avoid using variable names that begin with capital letters), I will name all matrices beginning with an m (so mA, mC, mWhatever) and all vectors beginning with a v (vX, vSolution, vSomething).

So the system becomes:

system = mA . vX == vC

Solve[system]

So we can solve the system very easily this way. 

Now, let's set up a system of differential equations.  Let's set up:
({{x' = x - 2  y}, {y' = 2  x + y}}) with initial conditions x(0) = 6,  y(0) = 4.

Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution]

mA=(({{1, -2}, {2, 1}}));vX[t_] = {x[t], y[t]}

system = vX'[t] == mA . vX[t]

vSolution[t_] = vX[t] /. DSolve[system, vX[t], t]

Or, if you prefer to include the initial conditions in the solution:

mX0 = {x0, y0}

vSolution[t_] = 
      vX[t] /. DSolve[{system, vX[0] == mX0}, vX[t], t]

A more complicated system might look something like this:

Clear[mA, vX, vC, x, y, t, x0, y0, vX0, vSolution]

mA=(({{1, -2}, {2, 1}}));vC = {2, t};vX[t_] = {x[t], y[t]}

system = vX'[t] == mA . vX[t] + vC

vSolution[t_] = vX[t] /. DSolve[system, vX[t], t]

Working with matrices and eigensystems

Sometimes, the solution that DSolve gives to a system isn't really what you were looking for (it may look much more complicated than you expected).  This is because Mathematica has to be able to solve some very complicated systems.  Sometimes, if your system is relatively simple (an autonomous homogeneous system, for example), you might be able to get a simpler solution by finding the eigenvalues and eigenvectors.  For example, if you are trying to solve the following system:

Clear[mA, vX, x1, x2, x3, x4, system, vSolution]

mA=(({{2, 3, -2, 1}, {5, 0, 6, -1}, {7, 5, 9, -7}, {8, 4, 2, 1}}));vX[t_] = {x1[t], x2[t], x3[t], x4[t]}

system = vX'[t] == mA . vX[t]

Let's find the eigenvalues and eigenvectors of this system (I stuck the N[...] in there to get a numeric result because the exact values are pretty insane...):

eigValues = Eigenvalues[N[mA]]

eigVectors = Eigenvectors[N[mA]]

Notice that, since all the coefficients of the original system are real, the complex ones occur in complex conjugate pairs.  Now, you can just cut and paste the appropriate eigenvalues/eigenvectors to construct the solution, or with a little more finesse (I'll take them from last to first; eigVectors[[4]] means the 4th vector in the list):

vX4 = 
      eigVectors[[4]] (E^(eigValues[[4]] t))

vX3 = 
      eigVectors[[3]] (E^(eigValues[[3]] t))

Now, the other two eigenvalues are complex, so recall that we handle those a bit differently:

vX2 = (E^(Re[eigValues[[1]]] 
              t)) (Re[eigVectors[[1]]] 
              Cos[Im[eigValues[[1]]] t] - 
            Im[eigVectors[[1]]] 
              Sin[Im[eigValues[[1]]] t])

vX1 = (E^(Re[eigValues[[1]]] 
              t)) (Re[eigVectors[[1]]] 
              Cos[Im[eigValues[[1]]] t] + 
            Im[eigVectors[[1]]] 
              Sin[Im[eigValues[[1]]] t])

(Wasn't that fun?  Make sure you understand where all of these came from before you continue.)  

Now, we can find the general solution by combining these:

vSolution[t_] = 
      C[1] vX1 + C[2] vX2 + C[3] vX3 + C[4] vX4 // FullSimplify

vSolution[t] // MatrixForm

If the initial conditions are given by:

vX0 = {1, 1, 1, 1}

We can solve for the constants:

initConditions = Solve[vSolution[0] == vX0]

The error here is a warning that we probably need to increase our numeric precision for all of our calculations (especially finding the eigenvalues and eigenvectors); round-off errors could easily screw things up for this complicated a situation.  However, we finally get:

vSolution[t] /. initConditions // Transpose // 
      MatrixForm

Of course, if you have eigenvalues with multiplicity greater than 1, then you have a lot more pain and suffering to deal with.  Unfortunately, Mathematica doesn't automatically give you all the eigenvalues ("generalized eigenvalues") you need for this, though you can certainly use Mathematica to find them.  However, there is a much better way to deal with systems of equations (with constant coefficients, anyway), which I cover in the next section, so let's just skip that, shall we?

Solving using exponentiation of matrices (or "exponentiation of operators")

When you solve the single differential equation x' = a x,  x(0) = (Subscript[x, 0]), you get the solution x(t) = (Subscript[x, 0]) (e^(a t)).  It turns out that if you are solving the system of differential equations X'=A .  X, X(0)=(Subscript[X, 0]) (where A is a matrix of constant coefficients), the solution winds up being X(t)=(e^(A t)).(Subscript[X, 0]) (needless to say, there is a lot of math that goes on behind the scenes to make this work out).  Of course, for this to make any sense, you have to have some idea what it means to raise e to a matrix power.  I won't go into all the details here, but the short form is that if you consider one definition of (e^a) to be (∑_(n = 0)^∞) ((a^n)/(n!)), then you can define (e^A) as (∑_(n = 0)^∞)((A^n)/(n!))(where (A^0)=I, the identity matrix), assuming the infinite sum actually converges.  So, for example, (e^((({{2, 0}, {0, 3}}))t))=(∑_(n = 0)^∞)(((((({{2, 0}, {0, 3}})))^n)t)/(n!))=(({{(∑_(n = 0)^∞)((((2  t))^n)/(n!)), 0}, {0, \(∑\+\(n = 0\)\%∞\((3  t)\)\^n\/\(n!\)\)}}))=(({{(e^(2  t)), 0}, {0, (e^(3  t))}}))  (this is, of course, a very simple example).  The good news is that Mathematica has a function built in for dealing with this: MatrixExp.  So, we could use:

MatrixExp[(({{2, 0}, {0, 3}}))t]//MatrixForm

So, let's work the problem from the previous section again, shall we?

Clear[mA, vX0, vSolution]

mA=(({{2, 3, -2, 1}, {5, 0, 6, -1}, {7, 5, 9, -7}, {8, 4, 2, 1}}));vX0 = {1, 1, 1, 1};

vSolution[t_] = MatrixExp[N[mA] t] . vX0 // ComplexExpand

(Notice that I use N[mA] in order to get a nicely readable solution, sort of...)

vSolution[t] // MatrixForm

Notice, however, the imaginary components in the answer.  This is due to some round-off errors.  Now, for positive t values, this turns out to be relatively insignificant.  Compare:

vSolution[100. ] // MatrixForm

However, for negative t values, it certainly isn't:

vSolution[-100. ] // MatrixForm

You can fix this in one of two ways:  The best way is to keep things symbolic as long as possible:

vSolution[t_] = MatrixExp[mA t] . vX0 // ComplexExpand

This isn't a very readable solution, but notice that:

vSolution[100. ] // MatrixForm

vSolution[-100. ] // MatrixForm

The other method is to just take the real components of your solution (and hope that nothing important is going on behind the scenes with that imaginary component):

vSolution[t_] = MatrixExp[N[mA] t] . vX0 // Re // 
        ComplexExpand

vSolution[100. ] // MatrixForm

vSolution[-100. ] // MatrixForm

General symbolic solutions

Of course, if you really don't want to fool with any of the above, you could just call up DSolve:

Clear[mA, vX0, vSolution]

mA=(({{2, 3, -2, 1}, {5, 0, 6, -1}, {7, 5, 9, -7}, {8, 4, 2, 1}}));vX0 = {1, 1, 1, 1};vX[t_] = {x1[t], x2[t], x3[t], x4[t]};

system = vX'[t] == mA . vX[t]

vSolution[t_] = 
      vX[t] /. First[
          DSolve[{system, vX[0] == vX0}, vX[t], t]]

Just for comparison:

vSolution[100. ] // MatrixForm

vSolution[-100. ] // MatrixForm

If you have a nonhomogeneous system:

mA=(({{2, 3}, {5, 0}}));vX0 = {1, 1};vX[t_] = {x1[t], x2[t]};vC = {Sin[t], Cos[t]};

system = vX'[t] == mA . vX[t] + vC

vSolution[t_] = 
      vX[t] /. First[
          DSolve[{system, vX[0] == vX0}, vX[t], t]]

Numeric approxiamtions

Or, we could have Mathematica do all the nasty approximations behind our backs and use NDSolve:

Clear[mA, vX0, vSolution]

mA=(({{2, 3, -2, 1}, {5, 0, 6, -1}, {7, 5, 9, -7}, {8, 4, 2, 1}}));vX0 = {1, 1, 1, 1};vX[t_] = {x1[t], x2[t], x3[t], x4[t]};

system = vX'[t] == mA . vX[t]

vSolution[t_] = 
      vX[t] /. First[
          NDSolve[{system, vX[0] == vX0}, 
            vX[t], {t, 0, 200}]]

This error message is a warning that you can't trust this approximation beyond about t = 21.  To fix that, you must increase the maximum number of steps (which will take longer to compute):

vSolution[t_] = 
      vX[t] /. First[
          NDSolve[{system, vX[0] == vX0}, vX[t], {t, 0, 200}, 
            MaxSteps -> 10000]]

Now, as before:

vSolution[100. ] // MatrixForm

vSolution[-100. ] // MatrixForm

You had better pay attention to this warning.  What is the problem?  We told NDSolve to solve on the domain 0 ≤ t ≤ 200, not for negative t.  This is something you need to watch out for.

Vector fields/direction fields/phase portraits

When you are working with systems of differential equations, it is often useful to look at either the vector field (direction field) of the system or the parametric plot of the solutions (or to just plot the individual solutions vs. time).  Let's take a look at how to do this with 2-dimensional systems (some of this is applicable to 3 dimensional systems as well).

(Sections with an * are optional.)

Vector fields

The vector field for a system of differential equations (first order, autonomous) is pretty easy to plot using PlotVectorField:

<< Graphics`

Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution]

mA=(({{1, -2}, {2, 1}}));vX[t_] = {x[t], y[t]}

system = vX'[t] == mA . vX[t]

PlotVectorField[mA . {x, y}, {x, -10, 10}, {y, -10, 10}, 
      Axes -> True]

(Notice, you don't have to actually solve the system to plot this.)  If you want to plot the "direction field", you just need to arrange that all the vectors have the same length:

directionField = 
      PlotVectorField[mA . {x, y}, {x, -10, 10}, {y, -10, 10}, 
        Axes -> True, ScaleFunction -> (0.4 &), 
        ScaleFactor -> 0.9]

Notice that these are plots in the XY plane (i.e., time doesn't enter into this).  So, what would it mean to plot solution curves on this plane?  It would mean to graph the parametric equations:  {x (t), y (t)}, graphing x vs. y in the plane (the "trajectory" in the XY plane).

vSolution[t_] = 
      vX[t] /. DSolve[{system, vX[0] == {x0, y0}}, vX[t], t] // 
        Flatten

If x(0) = 2 and y(0) = -1, then we have:

x0 = 2; y0 = -1

ParametricPlot[Evaluate[vSolution[t]], {t, 0, 2  π}]

To show a bunch of curves:

y0 = 1; Clear[x0]

solutionCurves = Table[vSolution[t], {x0, -5, 5, 1}]

solutionPlot = 
      ParametricPlot[Evaluate[solutionCurves], {t, -2 π, 2  π}, 
        PlotStyle -> Blue]

Show[directionField, solutionPlot, 
      PlotRange -> {{-10, 10}, {-10, 10}}]

Parametric plot vs. time plots

In interpreting the plots above, remember that these are parametric plots of X vs. Y.  What if you want to see how the two components (X and Y) behave over time?

Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution]

mA=(({{1, -2}, {2, 1}}));vX[t_] = {x[t], y[t]}

system = vX'[t] == mA . vX[t]

vSolution[t_] = 
      vX[t] /. DSolve[{system, vX[0] == {x0, y0}}, vX[t], t] // 
        Flatten

x0 = 2; y0 = -1

ParametricPlot[Evaluate[vSolution[t]], {t, 0, 5}]

xTimePlot = 
      Plot[Evaluate[vSolution[t][[1]]], {t, 0, 5}, 
        PlotStyle -> Red]

yTimePlot = 
      Plot[Evaluate[vSolution[t][[2]]], {t, 0, 5}, 
        PlotStyle -> Green]

Show[xTimePlot, yTimePlot]

(To get the first coordinate of vSolution[t], you use vSolution[t][[1]], and similarly for other coordinates.)How does this compare with the parametric plot shown earlier?

Nonhomogeneous systems and higher dimensional systems*

If your system is nonhomogeneous:

Clear[mA, vX, x, y, t, x0, y0, vX0, vSolution, vC]

mA=(({{1, -2}, {2, 1}}));vX[t_] = {x[t], y[t]};vC = {t, (t^2)};

system = vX'[t] == mA . vX[t] + vC

The problem is that this system doesn't just vary with X and Y, it also varies with t.  If you think of t as time (which it often is), however, you can use things like animation to show this.  One way to animate a graph is to make a table of "frames", which you can then double-click on and animate.  So, if we use one frame for every "tick" of the clock we can make some interesting vector fields (we let t go from 0 to 10, with a stepsize of 1 here):

Table[
      PlotVectorField[
        Evaluate[mA . {x, y} + vC], {x, -10, 10}, {y, -10, 10}, 
        Axes -> True, ScaleFunction -> (0.4 &), 
        ScaleFactor -> 0.9], {t, 0, 10, 1}]

While this is pretty cool, I'm not sure it is really all that useful.  (You can't, for example, use one of these graphs to trace out a trajectory.  You would have to trace it out in real-time while it is being animated...)

Notice, that we don't have to worry about this when graphing the parametric equations:  {x(t), y(t)}, graphing x vs. y in the plane (the "trajectory" in the XY plane).  It is just another parametric equation (though a bit ugly...)

vSolution[t_] = 
      vX[t] /. DSolve[{system, vX[0] == {x0, y0}}, vX[t], t] // 
        Flatten

If x(0) = 2 and y(0) = -1, then we have:

x0 = 2; y0 = -1;

ParametricPlot[Evaluate[vSolution[t]], {t, 0, 2  π}]

Notice how this has been "distorted" from the homogeneous case above.

So, what if you have a system of more than 2 dimensions?  Well, if it is three dimensional, you could use the 3-dimensional graphing capabilities of Mathematica to draw a vector field or trajectories.  (This may or may not be useful to you.)  Another possibility is to take the components 2 at a time and treat them like a 2-dimensional case (to tell you something at least, of how those 2 components interact).  I may put up some examples of this later if I have time.

Phase portraits of numeric approximations (think nonlinear systems)

Unfortunately, many nonlinear systems of differential equations can't be solved (by Mathematica, at least) in any reasonable sort of manner. You can use the critical points of the system (we are talking mainly about 2-dimensional systems here) along with the eigenvalues of the linear approximaiton to the system and its phase portrait to analyze these systems.  While the material given above is theoretically sufficient to draw these things, it turns out there are some annoying little issues that can really cuase problems if you aren't careful (and sometimes even if you are).

So, let's start with the following nonlinear autonomous system:

({{x' = (x^2) - 2  (y^2) - 2}, {y' = -4 (x^2) + 3  (y^2)}}) or ({{x' = f(x, y)}, {y' = g(x, y)}}) where f(x, y) = (x^2) - 2  (y^2) - 2 and g(x, y) = -4 (x^2) + 3  (y^2)

f[x_, y_] := (x^2) + 2  (y^2) - 2; 
    g[x_, y_] := -4 (x^2) + 3  (y^2)

Let's find our critical points:

critPts = 
      Solve[{f[x, y] == 0, g[x, y] == 0}] // N

Now, to classify these critical points, we need to evaluate the Jacobian at each of them and find its Eigenvalues (at the origin).  So, we need to define a Jacobian first:

jacob[f_, g_]:=(({{(∂_x) f, (∂_y) f}, {(∂_x) g, (∂_y) g}}))

jacob[f[x, y], g[x, y]]

At each critical point, you need to find the Jacobian and then find the eigenvalues of that matrix.  So, for example, at the first critical point:

jacob[f[x, y], g[x, y]] /. critPts[[1]] // MatrixForm

Eigenvalues[jacob[f[x, y], g[x, y]]] /. critPts[[1]]

So, telll me, what does it mean that these critical points are imaginary with negative real part?  (Chcck your answer down below when we graph the phase portrait.)  Now, it is actually pretty easy to do this for all of the critical points at once (each row is a pair of eigenvalues, listed in the same order as the critical points in the list above):

Eigenvalues[jacob[f[x, y], g[x, y]]] /. critPts // TableForm

Take a moment to go through and classify each critical point.  Now, let's graph the vector field (well, the direction field, to be precise) for this system:

<< Graphics`

vField = 
      PlotVectorField[{f[x, y], g[x, y]}, {x, -2, 2}, {y, -2, 2}, 
        Axes -> True, ScaleFunction -> (0.4 &), 
        ScaleFactor -> 0.2, 
        ColorFunction -> (Green  &)]

(Notice that we choose our domain to plot over so we can see the critical points and the region near them.  If your vectors are too long or too short, you can change the value of the ScaleFactor option.)

Let's graph our critical points:

critPtsPlot = 
      ListPlot[{x, y} /. critPts, 
        PlotStyle -> {Red, PointSize[0.02]}]

(Notice two things here:  We use ListPlot since we just want to plot some points, and we choose a PointSize of about 2% to make the points big enough to see easily.)

Sometimes, it is also handy to graph your "null-clines" as well:

nullClines = {f[x, y] == 0, g[x, y] == 0}

nullClinePlot = 
      ImplicitPlot[nullClines, {x, -2, 2}, {y, -2, 2}, 
        PlotStyle -> {Blue, Yellow}, PlotPoints -> 100]

Now, it's time to put them all together:

Show[vField, nullClinePlot, critPtsPlot]

This is pretty cool, but it would also be really helpful to add in some trajectories as well.  Unfortunately, there is really not a nice automatic way to pick out the most useful trajectories to plot.  Basically, you have to choose these individually (unless you just get lucky).  Basically, you pick points in interesting regions of the graph to be initial conditions.  Of course, since that means those would be the starting points of the graph, you might want to actually plot them for some negative t values as well as postive (i.e., backtrack along the curves to see where they came from).  

Now, if the system can be solved symbolically in a reasonable amount of time, this isn't all that difficult (see the examples earlier in this notebook).  However, if we try to solve this one, we get:

DSolve[{x'[t] == f[x[t], y[t]], y'[t] == 
          g[x[t], y[t]], x[0] == x0, y[0] == y0}, {x[t], y[t]}, 
      t]

So no dice there.  Clearly, we need to use NDSolve to find the solution.  So, for example, if we want to plot a curve with initial condition (1, 0) we could do:

ParametricPlot[
      Evaluate[{x[t], y[t]} /. 
          First[NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == 
                  g[x[t], y[t]], x[0] == 1, y[0] == 0}, {x[t], 
                y[t]}, {t, 0, 2}]]], {t, 0, 2}]

This sounds good.  But what about through the point (1, -1)?

ParametricPlot[
      Evaluate[{x[t], y[t]} /. 
          First[NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == 
                  g[x[t], y[t]], x[0] == 1, 
                y[0] == -1}, {x[t], y[t]}, {t, 0, 2}]]], {t, 0, 
        2}]

Uh, I think we have a problem here.  If you read the error messages, it is warning you that it can't compute a good solution all the way out to t = 2 (it actually doesn't get much beyond t = 0.5), so it tells you it is extrapolating to get all the way to 2.  As you can probably guess from looking at the numbers on the axes, this was pretty wildly unsuccessful.  You could fix this by simply graphing from t = 0 to t = 0.5:

ParametricPlot[
      Evaluate[{x[t], y[t]} /. 
          First[NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == 
                  g[x[t], y[t]], x[0] == 1, 
                y[0] == -1}, {x[t], y[t]}, {t, 0, 0.5}]]], {t, 0, 
        0.5}]

This works okay if you are only graphing a few curves, but it would really be nice to supply a list of points and then have Mathematica graph the solution curves for all of them, without having to set separate domains for individual curves.  The best way to do this would be to just instruct Mathematica to not extrapolate approximate functions beyond their domain.  Unfortunately, I haven't been able to figure out a way to do that (except on a case by case basis).  

So, after spending way too much time on this over the weekend, I decided to just write my own function to take care of this.  This was really annoying.  The good news for you is that you don't really have to worry about it.  You can just use the function I wrote (think cut-and-paste) and everything works out okay.  (You owe me for this one...)  Now, since I really didn't want to spend even more of my weekend making this really fancy, this means that you have to use this function exactly as I show you; I haven't built in any fancy error-handling or usage notes or anything else.  If you prefer not to use it, you can just plot individual curves (naming each one of them) and then combine them all together using the Show command.  On the other hand, if you prefer to take advantage of my function, I list it below (don't forget to execute it); do not change anything in it (unless you really know Mathematica pretty well) or you might break it.  (Unless you are just interested, it really isn't at all important that you understand how it works.  Actually, if you need to, you could change the value for PlotPoints, but be careful.)

TrajectoryPlot[system_, 
        initialPts_, {time_, tMin_, tMax_}, {xVar_, xMin_, xMax_}, {yVar_, 
          yMin_, yMax_}] := 
      Module[{}, 
Off[
          NDSolve::ndsz]; 
Off[
          InterpolatingFunction::dmval]; 

getSolution[k_] := 
          First[NDSolve[
              Join[system, {x[0] == initialPts[[k, 1]], 
                  y[0] == initialPts[[k, 2]]}], {x, y}, {t, tMin, 
                tMax}]]; 

getEstDomain[
            funct_] := {t, Flatten[funct[[1]]][[1]] + 
              0.2, Flatten[funct[[1]]][[2]] - 
              0.2}; 

Block[{$DisplayFunction = Identity}, 
solutionList = 
              Table[ParametricPlot[Evaluate[{x[t], y[t]} /. getSolution[k]], 
                  Evaluate[getEstDomain[x /. getSolution[k]]], 
                  PlotRange -> {{xMin, xMax}, {yMin, yMax}}, 
                  AspectRatio -> Automatic, PlotPoints -> 100], {k, 
                  1, Length[initialPts]}];]; 
On[
          NDSolve::ndsz]; 
On[
          InterpolatingFunction::dmval]; 

Show[solutionList]
]

So, how do you use it?  Fortunately, it is very simple.  You just have to supply it with the system (without any sort of initial conditions), a list of points you want curves through, the domain you want t to vary over, and the domains for x and y.  Like this:

points = {{1, 1}, {1, -1}, {-1, 
          1}, {-1, -1}}

solCurves = 
      TrajectoryPlot[{x'[t] == f[x[t], y[t]], y'[t] == 
            g[x[t], y[t]]}, 
        points, {t, -100, 100}, {x, -2, 2}, {y, -2, 
          2}]

Now, let's put all our graphs together:

Show[vField, nullClinePlot, critPtsPlot, solCurves, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Let's skip the null-clines and add in some more points (you should think about what points to choose to make things clearer):

points = {{1, 1}, {1, -1}, {-1, 
            1}, {-1, -1}, {0,  .25}, {-1, -2}, {0, -2}, {0,  .8}, {-1, 1.5}, {1.5, 1.5}, {- .5, 1}, { .5, -2}, {-1, 
            1.2}, {0.5, 0}};

solCurves = 
      TrajectoryPlot[{x'[t] == f[x[t], y[t]], y'[t] == 
            g[x[t], y[t]]}, 
        points, {t, -100, 100}, {x, -2, 2}, {y, -2, 
          2}]

Show[vField, critPtsPlot, solCurves, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Or, if you're just lazy (like me), you could have Mathematica plug in a grid of points for you (and if you're lucky, it will look good):

points = 
      Flatten[Table[{x, y}, {x, -2, 2, 0.5}, {y, -2, 2, 0.5}], 
        1]

(You need the Flatten command here to knock out some excess {}'s.)

solCurves = 
      TrajectoryPlot[{x'[t] == f[x[t], y[t]], y'[t] == 
            g[x[t], y[t]]}, 
        points, {t, -100, 100}, {x, -2, 2}, {y, -2, 
          2}]

Show[vField, critPtsPlot, solCurves, 
      PlotRange -> {{-2, 2}, {-2, 2}}]

Well, you get the idea.  The combination of the direction field with the solution curves actually gives you a pretty decent idea of what's going on here.  Of course, you still have to figure out what that means, but that depends on the specific application.