The procedure for solving a coupled system of differential equations follows closely that for solving a higher order differential equation. In fact, you can think of solving a higher order differential equation as just a special case of solving a system of differential equations.

**Systems of first order differential equations**

To solve a system of first order differential equations:

• Define a vector containing the initial values of each unknown function.

• Define a vector-valued function containing the first derivatives of each of the unknown functions.

• Decide which points you want to evaluate the solutions at

• Pass all this information into rkfixed.

The rkfixed function will return a matrix whose first column contains the points at which the solutions are evaluated and whose remaining columns contain the solution functions evaluated at the corresponding point, Figure 16-5 shows an example solving the equations:

**Systems of higher order differential equations**

The procedure for solving a system of nth order differential equations is similar to the procedure for solving a system of first order differential equations.

The main differences are:

• The vector of initial conditions must contain initial values for the n – 1 derivatives of each unknown function in addition to initial values for the functions themselves .

• The vector-valued function must contain expressions for the n – 1 derivatives of each unknown function in addition to the nth derivative.

The example in Figure 16-6 shows how to go about solving the system of second order differential equations

The function rkfixed returns a matrix in which:

• The first column contains the values at which the solutions and their derivatives are to be evaluated.

• The remaining columns contain the solutions and their derivatives evaluated at the corresponding point in the first column. The order in which the solution and its derivatives appear matches the order in which you put them into the vector of initial conditions that you passed into rkfixed.

**Specialized differential equation solvers**

The rkfixed function discussed thus far is a good general-purpose differential equation solver. Although it is not always the fastest method, the Runge-Kutta technique used by this function nearly always succeeds. Mathcad Professional includes several more specialized functions for solving differential equations, and there are cases in which you may want to use one of Mathcad’s more specialized differential equation solvers. These cases fall into three broad categories:

• Your system of differential equations may have certain properties which are best exploited by functions other than rkfixed. The system may be stiff (Stif.fb, Stiffr); the functions could be smooth (Bulstoer) or slowly varying (Rkadapt).

• You may have a boundary value rather than an initial value problem (sbval and bvalfit).

• You may be interested in evaluating the solution only at one point (bulstoer, rkadapt, stif.fb and stiffr).

You may also want to try several methods on the same differential equation to see which one works the best. Sometimes there are subtle differences between differential equations that make one method better than another.

The following sections describe the use of the various differential equation solvers and the circumstances in which they are likely to be useful

**Smooth systems**

When you know the solution is smooth, use the Bulstoer function instead of rkfixed. The Bulstoer function uses the Bulirsch-Stoer method rather than the Runge-Kutta method used by rkfixed. Under these circumstances, the solution will be slightly more accurate than that returned by rkfixed.

The argument list and the matrix returned by Bulstoer is identical to that for rkfixed

**Slowly varying solutions**

Given a fixed number of points, you can approximate a function more accurately if you evaluate it frequently wherever it’s changing fast and infrequently wherever it’s changing more slowly. If you know that the solution has this property, you may be better off using Rkadapt. Unlike rkfixed which evaluates a solution at equally spaced intervals, Rkadapt examines how fast the solution is changing and adapts its step size accordingly.

This “adaptive step size control” enables Rkadapt to focus on those parts of the integration domain where the function is rapidly changing rather than wasting time integrating a function where it isn’t changing all that rapidly.

Note that although Rkadapt will use nonuniform step sizes internally when it solves the differential equation, it will nevertheless return the solution at equally spaced points. Rkadapt takes the same arguments as rkfixed. The matrix returned by Rkadapt is identical in form to that returned by rkfixed.

**Evaluating only the final value**

The differential equation functions discussed so far presuppose that you’re interested in seeing the solution y(x) over a number of uniformly spaced x values in the integration interval bounded by xl and x2. There may be times, however, when all you want is the value of the solution at the endpoint, y(x2). Although the functions discussed so far will certainly give you y(x2), they also do a lot of unnecessary work returning intermediate values of y(x) in which you have no interest.

If you’re only interested in the value of y(x2) , use the functions listed below. Each function corresponds to one of those already discussed. The properties of each of these functions are identical to those of the corresponding function in the previous sections

**Boundary value problems**

So far, all the functions discussed in this chapter assume that you know the value taken by the solutions and their derivatives at the beginning of the interval of integration. In other words, these functions are useful for solving initial value problems.

In many cases, however, you may know the value taken by the solution at the endpoints of the interval of integration. A good example is a stretched string constrained at both ends. Problems such as this are referred to as boundary value problems. The first section

discusses two-point boundary value problems: one-dimensional systems of differential equations in which the solution is a function of a single variable and the value of the solution is known at two points. The section following this discusses the more general case involving partial differential equations.