Contents Previous Next Subchapters Current Chapters-> trapz gaussq quad2d gaussq2d gaussleg quadint int2inf ode4rk expmat ode45rk odepade odestiff Parent Chapters-> Omatrix6 integration odestiff Search Tools-> contents reference index search

Solving Stiff Ordinary Differential Equations
 Syntax `odestiff(function `fval`, `order`, `ti`, `tout```, ...          ```maxstep`, `minstep`, `yi`, `eps`, `skip`)` See Also ode4rk , odepade

Description
Uses a backward difference method to solve the differential equation ```      d      -- y(t) = f(t, y(t))      dt ```with initial condition `y(t) = yi` at `t = ti`. The integer scalar order specifies the order of the backward difference method. The real or double-precision scalar ti specifies the initial value for `t`. The row vector tout has the same type as ti and specifies the values of `t` at which `odestiff` should calculate the value of `y`. The scalar maxstep has the same type as ti and specifies the maximum step size in `t`. The scalar minstep has the same type as ti and specifies the minimum step size in `t`. The column vector yi has the same type as ti and specifies the initial value of `y`. The column vector eps has the same type and dimension as yi and `eps(i)` is the desired accuracy for the i-th element of `y(t)`. The integer scalar skip specifies the level of tracing inside of `odestiff`. If `skip = 0`, no tracing is done; otherwise, the value of `y(t)` is printed at time `tout(skip)`, `tout(2 skip)` ... . The return value for `odestiff` has the same type and row dimension as yi, the same number of columns as tout, and its j-th column is the value of `y` at time `tout(j)`. ``` ```fval`(`tin`, `yin`, `fout`, `dfout```) ```This function call sets fout the value of `f(tin, yin)` and dfout to the value ```      d      -- f(tin, yin)      dy ```The scalar tin has the same type as ti. The vector yin has the same type and dimension as yi. The input value of fout does not matter. Its output value has the same type and dimension and row dimension as yi. The input value of dfout does not matter. Its output value is a square matrix with the same type and row dimension as yi and its (i,j)-th element is the partial of the i-th element of `f(t, y)` with respect to the j-th element of `y`.

Example
The following program plots the solution of the differential equation ```      y'(t) =  y (t)     y (0) = 0       1        2         1            y'(t) = -y (t)     y (0) = 1       2        1         2 ```The solution is ```      y (t) = sin(t)       1      y (t) = cos(t)       2 `````` clear # function fval(tin, yin, fout, dfout) begin      fout   = {yin(2), -yin(1)}      dfout  = {[0., 1.], [-1., 0.]}      return end order   = 4 ti      = 0. tf      = 2 * pi maxstep = 1. minstep = .01 yi      = {0., 1.} eps     = {1e-4, 1e-4} skip    = 0 tout    = seq(30)' * (tf - ti) / real(30) yout    = odestiff(function fval, order, ti, tout, ...                    maxstep, minstep, yi, eps, skip) gplot(tout', yout') ```