Nothing

```
odesolv <- function(bwtlist, ystart=diag(rep(1,norder)),
h0=width/100, hmin=width*1e-10, hmax=width*0.5,
EPS=1e-4, MAXSTP=1000)
{
# Solve L u = 0,
# L being an order M homogeneous linear differential operator,
# (Lu)(t) = w_1(t) u(t) + w_2(t) Du(t) + ...
# w_m(t) D^{m-1} u(t) + D^m u(t) = 0
# for function u and its derivatives up to order m - 1.
# Each such solution is determined by the values of u and its
# m - 1 derivatives at time t = 0. These initial conditions are
# contained in the columns of the matrix YSTART, which has exactly
# m rows. The number of solutions computed by ODESOLV is equal to
# the number of columnsof YSTART. In order for the solutions to be
# linearly independent functions, the columns of YSTART must be
# linearly independent. This means that the maximum number of
# linearly independent solutions is m.
# The solution for each value of t is a matrix, y(t).
# Any column of y(t) contains u, Du, ... , D^{m-1} u at argument
# value t, for the corresponding set of starting values for these
# first m derivatives. It is the job of this function to estimate
# these values, and ODESOLV will choose a set of values of TNOW at
# which these can be estimated with satisfactory accuracy.
# ODESOLV uses the Runge-Kutta method, which is a good general
# purpose numerical method. But it does not work well for stiff
# systems, and it can fail for poor choices of initial conditions
# as well as other problems.
# Arguments:
# BWTLIST ... a list containing m = 1 weightfunctions.
# The weight functions w_1, ... , w_m are the functions with
# indices 2, 3, ..., m+1.
# YSTART ... initial values for Y. This is a matrix with M rows,
# were M is the order of the operator L. Any column of M
# specifies intial values for derivatives 0, 1, ... M-1. Each
# column must specify a unique set of initial conditions. A
# frequent choice is the identity matrix of order M, and this is
# the default.
# H0 ... initial trial step size
# HMIN ... minimum step size
# HMAX ... maximum step size
# EPS ... error tolerance
# MAXSTP ... maximum number of Runge Kutta steps permitted. If the equation
# is difficult to solve, this may have to be increased.
# Returns:
# TP ... vector of TRUE values used
# YP ... m by m by length(TP) array of Y-values generated for
# values in TP.
# Note that ODESOLV calls function DERIVS in order to evaluate the
# differential operator. Also, it works by redefining the order m
# linear differential equation as a linear system of m first order
# differential equations. DERIVS evaluates the right side of this
# linear system.
# Last modified 30 October 2005
MAXWARN <- 10
# determine the order of the system m
norder <- length(bwtlist)
# determine the range of values over which the equation is solved
bfdPar1 <- bwtlist[[1]]
wbasis <- bfdPar1$fd$basis
rangeval <- wbasis$rangeval
tbeg <- rangeval[1]
tend <- rangeval[2]
width <- tend - tbeg
tnow <- tbeg
h <- min(c(h0, hmax))
tp <- tnow
# set up the starting values
ystartd <- dim(ystart)
if (ystartd[1] != norder) stop("YSTART has incorrect dimensions")
n <- ystartd[2]
yp <- c(ystart)
index <- abs(yp) > 1e-10
y <- ystart
# Since ODESOLVE is slow, progress is displayed
cat("Progress: each dot is 10% of interval\n")
# initialize the solution
tdel <- (tend - tbeg)/10
tdot <- tdel
iwarn <- 0
# solve the equation using a maximum of MAXSTP steps
for (nstp in 1:MAXSTP) {
# evaluate the right side at the current value
dydt <- derivs(tnow, y, bwtlist)
yscal <- c(abs(y) + abs(h*dydt) + 1e-30)[index]
if (nstp > 1) {
tp <- c(tp,tnow)
yp <- c(yp,c(y))
}
if (tnow >= tdot) {
cat(".")
tdot <- tdot + tdel
}
if ((tnow+h-tend)*(tnow+h-tbeg) > 0) h <- tend-tnow
# take a Runge-Kutta step to the next value
result <- rkqs(y, dydt, tnow, h, bwtlist, yscal, index, EPS)
tnow <- result[[1]]
y <- result[[2]]
h <- result[[3]]
hnext <- result[[4]]
# test to see if interval has been covered, and return
# if it has.
if ((tnow-tend)*(tend-tbeg) >= 0) {
cat(".")
tp <- c(tp,tnow)
yp <- c(yp,c(y))
yp <- array(yp,c(norder,n,length(tp)))
return(list(tp, yp))
}
# test if the step is too small.
if (abs(hnext) < hmin) {
warning("Stepsize smaller than minimum")
hnext <- hmin
iwarn <- iwarn + 1
if (iwarn >= MAXWARN) stop("Too many warnings.")
}
# test if the step is too large.
h <- min(c(hnext, hmax))
}
warning("Too many steps.")
}
# ---------------------------------------------------------------
rkqs <- function(y, dydt, tnow, htry, bwtlist, yscal, index, EPS)
{
# Take a single step using the Runge-Kutta procedure to
# control for accuracy. The function returns the new
# value of t, the value of the solution at t, the step
# size, and a proposal for the next step size.
h <- htry
# check the accuracy of the step
result <- rkck(y, dydt, tnow, h, bwtlist)
ytemp <- result[[1]]
yerr <- c(result[[2]])[index]
errmax <- max(abs(yerr/yscal))/EPS
# modify the step size if ERRMAX is too large
while (errmax > 1) {
htemp <- 0.9*h*(errmax^(-0.25))
h <- max(c(abs(htemp),0.1*abs(h)))
tnew <- tnow + h
if (tnew == tnow) stop("stepsize underflow in rkqs")
result <- rkck(y, dydt, tnow, h, bwtlist)
ytemp <- result[[1]]
yerr <- result[[2]]
errmax <- max(abs(yerr/yscal))/EPS
}
# set up the proposed next step
if (errmax > 1.89e-4) {
hnext <- 0.9*h*(errmax^(-0.2))
} else {
hnext <- 5.*h
}
tnow <- tnow + h
y <- ytemp
return( list(tnow, y, h, hnext) )
}
# ------------------------------------------------------------------
rkck <- function(y, dydt, tnow, h, bwtlist)
{
# Take a single Runge-Kutta step.
# Return the solution, and an estimate of its error.
C1 <- 37/378
C3 <- 250/621
C4 <- 125/594
C6 <- 512/1771
DC5 <- -277/14336
DC1 <- C1 - 2825/27648
DC3 <- C3 - 18575/48384
DC4 <- C4 - 13525/55296
DC6 <- C6 - 0.25
ytemp <- y + h*0.2*dydt
ak2 <- derivs(tnow+0.2*h, ytemp, bwtlist)
ytemp <- y + h*(0.075*dydt + 0.225*ak2)
ak3 <- derivs(tnow+0.3*h, ytemp, bwtlist)
ytemp <- y + h*(0.3*dydt - 0.9*ak2 + 1.2*ak3)
ak4 <- derivs(tnow+0.6*h, ytemp, bwtlist)
ytemp <- y + h*(-11*dydt/54 + 2.5*ak2 + (-70*ak3+35*ak4)/27)
ak5 <- derivs(tnow+h, ytemp, bwtlist)
ytemp <- y + h*(1631*dydt/55296 + 575*ak2/512 + 575*ak3/13824 +
44275*ak4/110592 + 253*ak5/4096)
ak6 <- derivs(tnow+0.875*h, ytemp, bwtlist)
yout <- y + h*(C1*dydt + C3*ak3 + C4*ak4 + C6*ak6)
yerr <- h*(DC1*dydt + DC3*ak3 + DC4*ak4 + DC5*ak5 + DC6*ak6)
return (list( yout, yerr ) )
}
```

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.