Figure 1: Demonstrating the Levitron™
in suspension above the magnetic base.
The Mathematics
Summary
of Mathematical Foundations
Construction
of the Equations
Energy Equations
Appendix A- Driver Program for Runge-Kutta Solver
We are team 72 from Sandia Preparatory School. The
goal of our project has been to model the behavior of the Levitron™, a
small magnetic top, marketed as a child’s toy. The Levitron™ is capable
of hovering as it spins above a larger magnetic base. The interaction of
the magnetic fields of the top and the base, as well as the top’s rotational
motion, allow for it to sustain levitation for up to two minutes. We chose
this seemingly simple toy as the subject for our project because when one
delves into the mathematics and physics behind its behavior, one finds
that it actually exhibits chaotic properties. It has also been discovered
that subatomic particles trapped in magnetic fields demonstrate behavior
similar to that of the Levitron™. We chose to simulate its behavior because
our team had a mutual interest in chaos theory and the Levitron™ would
give us an introduction to some of the basics behind this newly emerging
field such as bifurcation and strange attractors. Chaos theory has become
a significant field because it can be applied to virtually any dynamic
system including weather patterns or the stock market. In an increasingly
complicated world, chaos theory has redefined this complexity and has the
potential to become one of the most significant foundations for future
scientific studies.
As most people who have played with magnets
before know, it seems to be impossible for one magnet to levitate above
another. The Levitron™ accomplishes this feat by virtue of its rotational
motion. The rotation of the top creates a gyroscopic motion that keeps
it floating as long as there is enough rotational momentum.
We are modeling its behavior mathematically through
Hamilton’s equations of motion, tracking the changes that occur over a
given time step. As explained in the body of this paper, we solve 12 differential
equations simultaneously using a method known as Runge-Kutta to accommodate
the time intervals between steps. Our code outputs sets of values for a
single variable. By graphing these variables with respect to time, we have
been able to study changes along the ‘x’ and ‘y’ axes and determined that
they act as a ‘strange attractor’. We have also modeled the effects of
momentum and are currently in the process of examining the chaotic behavior
of the angles of rotation and how they relate to each other. The top’s
behavior has proved to be an interesting subject that has potential significance
for future research on chaos.
Figure 1: Demonstrating the Levitron™
in suspension above the magnetic base.
It is important to discuss dynamic systems before
delving into the chaotic behavior that many demonstrate. A dynamic system
is any system that changes over time. This change can have a distinct pattern
or can be completely random. Virtually anything around us that involves
change can be classified as a dynamic system. Examples can include everything
from weather patterns to traffic patterns, or even the currents in the
ocean. Chaotic and random systems emerge mathematically because of their
extreme sensitivity to initial conditions. This separates them from more
ordered systems; a common example being a Newtonian physical problem in
which the answers do not change every time you solve the equations. We
can study chaotic systems through graphs that represent the solutions to
the system as the initial conditions evolve. Such graphs can include plots
of average traffic speed with respect to time, amount of moisture with
respect to time, or wind speed with respect to time. Creating such plots
allows the system to be studied graphically and lets us view
its patterns. Even such influential and everyday things as the stock market
demonstrate chaotic behavior.
A chaotic system is a type of dynamic system that
can be predicted by finding patterns in the solutions to its equations,
although a specific solution never repeats. In order to study these systems,
we solve several equations over time given certain initial conditions and
then plot the results on graphs. The chaotic patterns emerge on these graphs
after solving the equations with multiple sets of initial conditions. Changing
the conditions helps in the study of the behavior of these systems by allowing
us to perceive how each condition affects the system as a whole. The Levitron™
is a good example of this. Factors such as spin rate, height and position
along the ‘x’ and ‘y’ axes greatly effect both its stability and behavior.
There are different types of chaos. “Strange attractors”
are one type. When looking at a graph of a chaotic system, for example:
precession with respect to time, it will appear as though the graph is
trying to get to a certain point but never quite makes it. Strange attractors
often result in graphs that look like spirals but without perfect geometric
qualities. Bifurcation is another type of chaos that can be seen in the
Levitron™. Bifurcation is essentially the multiplication of frequencies,
which results in a branching pattern that perpetually increases in complexity.
Bifurcation is evident in the Levitron™ upon examining the frequency of
the ‘nodding’ or ‘wobble’ of the top itself. This is usually referred to
as nutation. The nutation may initially begin with a simple frequency and
then appear to become completely random. In fact, additional frequencies
are emerging and basically stacking upon one another, making the motion
increasingly complex as the top precesses about the spin axis. If more
than one nutation frequency is present in our model, it starts to divide
exponentially becoming four frequencies, then eight, and so on.
The Levitron™ is a spinning magnetic top that floats
above a square base with a ring magnet inside. Certain very precise conditions
must be met in order for the system to be stable enough to maintain prolonged
levitation. The magnetic field of the base has to be aligned parallel with
gravity or else the top will be forced to one side or the other. The rate
of spin of the top can neither be too slow nor too fast. If it is spun
too slowly, the torque caused from gravity cannot be overcome and the top
collapses. If it is spun too quickly, the top maintains an axis of rotation
which is too close to the vertical (in other words the tilt is not great
enough) and stability cannot be achieved. This tilt is known as precession
and is a crucial factor in our model. Without it, the top is prone to collapse.
Once stability is achieved, the top will float for a maximum of two minutes
because air resistance slows its rotation. The top can sustain levitation
for much longer periods in a vacuum. Some vacuum experiments have been
carried out, demonstrating that the top will spin for over 30 minutes.
The top eventually collapses because the multiplying nutation frequencies
increase the angle of precession leading to instability (this will be explained
further below). The only way to counteract this is to provide an outside
energy source that constantly re-aligns the top to the center of the field.
One such device is marketed as the ‘Perpetuator.’ This device keeps the
top rotating at the ideal spin rate of about 30 revolutions a second
and prevents the precession angle from becoming too large by the utilization
of another rotating magnet. The ‘Perpetuator’ can maintain stable levitation
for approximately four hours. The eventual collapse is usually because
of changes in the room temperature, which the Levitron™ is sensitive to
(this will further explored below as well).
The height above the base in which the top is stable
is very small: about 1.25 to 1.75 inches . Meeting these conditions is
accomplished by, first, spinning the top as close to the center of the
base as possible so that it can align itself in the center of the magnetic
field. Secondly, the lifter plate is brought up slowly until it is at the
proper height, at which point the top “jumps” off the plate and stabilizes
itself. The weight of the top also has to be very precise in order to balance
out the force of gravity and magnetic repulsion. Several thin, metallic
washers are used to change the weight of the top, which must be adjusted
to within a tenth of a gram. Slight temperature changes to the magnetic
base (plus or minus six degrees Fahrenheit) affect the magnet field enough
that the weight must be adjusted, thus excessive handling or even changes
in room temperature can effect its behavior .
You have seen that the stability of the top is dependent on several
factors, however we have not yet examined why these conditions are so sensitive.
From Earnshaw’s theorem, we discover that a stationary magnet cannot levitate
above another stationary magnet. One property of all bodies is that they
prefer to exist in the lowest possible energy state. For stability to be
achieved there must be a suitable minimum energy state.
The Mathematics:
Section
3.1 “A Summary of the Mathematical Foundations”:
The mathematics involved in the examination of the
conditions that govern the stability and motion of the Levitron™ can be
described through a set of equations that we shall attempt to describe
in the following section. Beginning with certain initial conditions and
assumptions, we can mathematically create a stable environment for the
spinning top and proceed to solve for the physical characteristics of the
system. Essentially, we begin with conditions for height, strength of the
dipole, and other such attributes that we know will be conducive to allowing
the top to levitate. This allows us to avoid having to 're-discover' exactly
where to begin since the region of stability for the top is so narrow.
Beginning with these conditions, we can move on
to examine the equations which stem from the standard equations for the
kinetic and potential energy of a spinning body. Using these equations
and a method known as nondimensionalization, we construct what is commonly
called the Hamiltonian equation. Nondimensionalization allows us to simplify
the equations by eliminating certain physical characteristics of the system
and treating the top as a single dipole without mass (this will be further
explained below). The nondimensionalized form of the Hamiltonian can then,
in turn, be used to establish a set of equations that solve for the necessary
variables which will be described in more depth. This set of motion equations
is analyzed over a given time interval to give values for the variables.
The initial conditions are then altered and the process is repeated. The
changes to the solutions can be graphed so that we can analyze them visually.
It is through this analysis that the chaotic behavior of the Levitron™
becomes evident.
Section 3.2 "Construction
of the Equations":
The first important mathematical concept in our
problem is the understanding of the two vectors that represent the position
and momentum of the system. Position (the ‘q’ vector) contains six variables
as shown below.
represents the arc that the top covers as it rotates and 
represents the motion of the top itself around its individual axis (see
figure 3.1).
Section 3.3 "Energy Equations”:
Once we have established the two primary vectors
in our problem, we can begin to construct the essential equations that
will allow us to evaluate these variables.
We begin with the equations for kinetic and potential
energy for a rotating body, which we are familiar with from basic physics.
We know that kinetic energy is equal to 1/2 the mass times the velocity
squared (1/2mv^2) in the simplest mechanical context. In this case, we
expand this to contain the necessary attributes of the spinning top:
’
contain the dimensional variables, including mass. The values chosen for
the constants are based on scientific observations of the Levitron™’s behavior.
Nondimensionalization is one way to arrange the essential elements while
eliminating a clutter of unnecessary variables and representing them as
more efficient constants. This not to say that the results are any less
accurate, we are merely utilizing this particular numerical method for
simplicity’s sake.
In order to solve the twelve equations presented
above, we will be utilizing the Runge-Kutta method. Runge-Kutta works by
transforming the differential equation in question to its slope form. Both
sides of the equation are multiplied by the change in ‘X’ or, in our case,
the change in ‘T’ (time). This allows us to solve for the current variable
by plugging in a step size for the change in ‘T’ and a previous value for
the variable we are solving for. This process is repeated to obtain sets
of values for the variable.
This is the simplest form of the Runge-Kutta method
and is known as a ‘First-Order Runge-Kutta Method’ or ‘Euler’s Method’.
This particular method has many limitations including a tendency to become
unstable and a relatively high error. To counter this, Runge-Kutta methods
have been developed that use more than one previous value to predict the
behavior of the equation in the future. This approach is called a ‘Higher-Order
Runge-Kutta Method’. This has the benefit of being relatively stable and
fairly accurate. The problem is that for each previous value that we consult
to get the current one, we add computation time. With the relatively small
time step we are using for our problem, this computation time adds up quickly;
this is why we need the super computers.
Runge-Kutta relies heavily on initial conditions.
For the first time step we must define each variable. As we have said before,
the goal of our project was not to study when the Levitron™ is stable but
to study its chaotic behavior. We want initial conditions that provide
for stable and continuos levitation. Several papers have been written about
the Levitron™ already, focusing mainly on stability. We can use their results
for the initial conditions and constants to create an environment suitable
for stable levitation. This allows us to build on previous results and
study the chaotic behavior -- an area left relatively unexplored. The keys
to stable levitation are the height of the top and the momentum of the
rotation. The Levitron™ is only stable at one height and with a high enough
momentum to counteract the torque from the base magnet. If we do not meet
these conditions, the system will collapse. Another interesting factor
is that if the momentum is too high, the system will also destabilize;
since the user controls this entirely, it is difficult to achieve just
the proper speed and twist when releasing the top. This only reiterates
how difficult it is to create stable levitation if the physical parameters
are not met. The height of the top, q3, is initialized at 1.72 inches.
The tilt, q4, is small (usually initialized between 0.001 and 0.015 radians)
. The momentum of the top itself, p6, is set between 5.0 and 9.0 radians
per second . This momentum is constant. The last value we must initialize
is the momentum around the vertical axis, p5, which is usually between
0.5 and 3.5 radians per second . We fixed the physical constants at a=0.089,
c=0.139, and 
=8.20; which
are consistent with observed values for the Levitron™ . The final facet
to be decided, after the initial conditions and constant values have been
dealt with, is the time step. With the nondimensionalization of the equation,
the consistent time scale is equal to which is the time it
would take the top to fall a distance of R/2 (for our model R=34.7mm so
time is equal to 59.5ms.) One thousand units of this time scale are equal
to approximately one minute.
One final note about our model: We do not take into
account the effects of air resistance. The only thing that air resistance
does is to slow down the rate of spin leading to collapse of the system.
We can reasonably assume that air resistance does not change the chaotic
behavior of the top in any great degree. Since our goal was to study the
chaotic behavior, we did not see the necessity of adding unneeded complexity.
We used several different Runge-Kutta routines,
as well as different computers. The first method was in ‘Waterloo Maple
V’. We programmed the equations into Maple V in order to solve small
problems on a PC that could be used to verify the results we receive from
our actual C code. The Maple V code functioned well but we were
unable to solve anything but the simplest problems due to the PC’s inability
to process a low enough error tolerance. The second option we used was
from the book, Numerical Recipes in C. We have programmed a driver that
utilizes the 5th order Runge-Kutta solver available in this book. The results
obtained from the program are compared to the simple problem we solved
using Maple V, and several test cases were compared to previously
reported results.
We are graphing the data using ‘Dataplot’ from Fortner
software. We have included three of these plots. The first two demonstrate
how a low initial momentum causes instability. Figure 4.1 shows a graph
of tilt versus time units. Remember that our time unit is 59.5ms so the
graph represents about half a second of real time. Initially the precession
of the top is within a very small range (between 0.001 to 0.015 radians)
but the spin rate is too low to counter the torque from gravity, causing
the tilt angle to increase at the 8th time unit. The increased tilt angle
is too high for the Levitron™ to continue levitating and the top collapses.
Figure 4.2 shows a graph of height versus time units and you can see how
the collapse of the top matches up with the increase in the tilt angle.
The last graph represents part of the top’s chaotic behavior. Figure 4.3
is a graph of the changes along the ‘x’ and ‘y’ axes. The chaotic behavior
that this position graph demonstrates shows one of the types mentioned
in the ‘Dynamics section,’ a strange attractor. If you look carefully you
can see the spiral-like change in the solutions. The graph represents 165
time units (approximately 10 seconds). We initialized the momentum of the
top, p6, to 8.7 radians a second; this is very close to the maximum value
under stable levitation. We are currently studying the chaotic behavior
of the angles and how they relate to time as well as to each other. These
relationships are where the truly interesting chaos emerges.
Figure 4.1: This graph shows tilt versus time
Figure 4.2: This shows height versus time.
Figure 4.3: This shows the changes in the ‘x’
and ‘y’ position.
Conclusions:
We have successfully modeled the Levitron™
for both unstable and stable configurations. We verified are results with
previously reported results on the stability of the Levitron™ and by comparing
the results to simpler solutions completed using Maple V on a PC.
Our most original accomplishment this year has been
the study of the chaotic behavior of the Levitron™. As we said earlier
in the paper, studies have been conducted on the Levitron™ in the past
but have focused almost exclusively on stability. We have primarily looked
at the chaotic behavior, building upon most previous work.
What remains to be done includes further study on
the chaotic behavior between the angles of rotation. We are currently attempting
to model the Levitron™ during sustained levitation in order to graph the
changes in the angles. The chaotic behavior emerges in these systems when
they are taken with respect to time and with respect to each other.
We have learned a great deal over the course of
the project, including how to work with non-linear mathematics. In addition,
we have gained a lot of practice in the area of computing with initial
value problems. We have come to understand the basics of chaotic theory
as well as the value in studying this field.
We would like to thank everyone at Los Alamos National
Labs for providing us the resources to complete these projects and for
their support over the course of our high school careers. It has been a
very positive program.
We owe the most gratitude to our advisor of four
years, Dr. Kevin Malloy, who has given up countless hours to help us understand
the math and science behind our projects. We cannot forget his colleagues
at the Center for High-Technology Materials at the University of New Mexico
who have provided technical assistance as well as access to software packages.
We extend a special thanks to Dr. Fred Gelbard from Sandia National Laboratories
who helped us debug our program and who was there to answer questions when
our advisor was not available. Last but certainly not least, we wish to
thank our sponsor, Mr. Neil McBeth who originally introduced us to the
program four years ago and who has always been there for us.
Berry, Michael V. The Physics of the Levitron. (Bibliographic
data unknown. From the
pamphlet that is included
with the Levitron.)
Chua, Leon; Madan, Rabinder N.. Sights and Sounds of Chaos. New
York: IEEE Service
Center, 1988. (Originally
published in IEEE Circuits and Devices Magazine).
Dullin, Holger R.; Easton, Robert W.. Stability of the Levitron.
Boulder, Colorado:
University of Colorado,
April 23, 1998.
Gans, Roger; Jones, Thomas B.; Washizu, Masao. Dynamics of the Levitron.
UK: IOP
Publishing Ltd., 1998.
Press, William; Flannery, Brian; Teukolsky, Saul; Vetterling, William.
Numerical
Recipes
in C: The Art of Scientific
Computing. Cambridge: Cambridge University Press, 1988.
Simon, Martin D.; Heflinger, Lee O.; Ridgeway, S.L..Spin Stablized
Magnetic Levitation.
California: UCLA Dept. of Physics (published in
the American Journal of Physics), unknown pub. date.
As well as the following software packages:
Waterloo Maple V
Fortner Software’s “DataPlot”
Driver Program for Runge-Kutta Solver
#include <stdio.h>
#include "nr.h"
#include "nrutil.h"
#include <math.h>
#define N 12 /*The number of simultaneous equations to solve*/
#define ITER 100
#define a 0.089
#define c 0.139
#define M 8.1880
#define MAXSTP 10000
#define TINY 1.0e-30
int kmax=0,kount=0; /* defining declaration */
float *xp=0,**yp=0,dxsav=0; /* defining declaration */
void derivs(x,y,dydx) /*Function to solve the right hand side of the
equations
used by odeint. DyDx is an array that holds the answers for each variable
represented the y array.*/
float x,y[],dydx[];
{
dydx[1] = y[7];
dydx[2] = y[8];
dydx[3] = y[9];
dydx[4] = (y[10]/a);
dydx[5] =0
(y[11]-(y[12]*cos(y[4])))/(pow(sin(y[4]),2.0));
dydx[6] =
(y[12]/c)-(((y[11]-(y[12]*cos(y[4])))*((cos(y[4])/sin(y[4]))*(1/sin(y[4]))))/a);
dydx[7] =
(1/(4*pow(((1.0+pow(y[3],2.0))),(9.0/2.0))))*(3*M*(-y[1]*cos(y[4])*(3-24*pow(y[3],2.0)
+8*pow(y[3],4.0))+y[3]*cos(y[5])*sin(y[4])*(-3-pow(y[3],2.0)+2*pow(y[3],4.0))));
dydx[8]=
(1/(4*pow((1.0+(pow(y[3],2.0))),(4.5))))*(3*M*(-y[2]*cos(y[4])*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0))
+y[3]*sin(y[5])*sin(y[4])*(-3-pow(y[3],2.0)+2*pow(y[3],4.0))));
dydx[9]=
(1/(8*pow((1.0+(pow(y[3],2.0))),(5.5))))*(-3*M*y[3]*cos(y[4])
*(-12-75*pow(y[2],2.0)-16*pow(y[3],2.0)+200*pow(y[2],2.0)*pow(y[3],2.0)
+4*pow(y[3],4.0)-40*pow(y[2],2.0)*pow(y[3],4.0)+8*pow(y[3],6.0)-5*pow(y[1],2.0)
*(15-40*pow(y[3],2.0)+8*pow(y[3],4.0))+(2*(1.0+pow(y[3],2.0))
*(2*pow((1.0+pow(y[3],2.0)),(.5))+8*pow(y[3],2.0)*pow((1.0+pow(y[3],2.0)),(.5))
+12*pow(y[3],4.0)*pow((1.0+pow(y[3],2.0)),(.5))+8*pow(y[3],6.0)*
pow((1.0+pow(y[3],2.0)),(0.5))+2*pow(y[3],8.0)*pow((1.0+pow(y[3],2.0)),(0.5))
-3*M*y[1]*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0))*cos(y[5])*sin(y[4])
-9*M*y[2]*sin(y[4])*cos(y[5])+2*pow(y[3],8.0)*pow((1.0+pow(y[3],2.0)),(0.5))
-3*M*y[1]*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0))*cos(y[5])*sin(y[4])
-9*M*y[2]*sin(y[4])*sin(y[5])+72*M*y[2]*pow(y[3],2.0)*sin(y[4])*sin(y[5])
-24*M*y[2]*pow(y[3],4.0)*sin(y[4])*sin(y[5])))));
dydx[10]=
(1/2)*(((2*y[12]*(y[11]-y[12]*cos(y[4]))/a*sin(y[4]))*((2*pow((y[11]
-y[12]*cos(y[4]),2.0)*cos(y[4]))/a*pow(sin(y[4]),2.0)*sin(y[4]))
-(1/(4*pow((1.0+(pow(y[3],2.0))),(4.5))))*(M*((-4-9*pow(y[2],2.0)+72*pow(y[2],2.0)
*pow(y[3],4.0)+8*pow(y[3],6.0)-3*pow(y[1],2.0)*(3-24*pow(y[3],2.0)+8*pow(y[3],4.0)))
*sin(y[4])-6*y[3]*(-3-pow(y[3],2.0)+2*pow(y[3],4.0)*cos(y[4])*(y[1]*cos(y[5])
+y[1]*sin(y[5]))))))));
dydx[11]=
(3*M*y[3]*(-3+2*pow(y[3],2.0))*sin(y[4]))*(y[2]*cos(y[5])-y[1]*sin(y[5]))/
((4*pow((1.0+(pow(y[3],2.0))),(3.5))));
dydx[12]=0;
}
void odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqc) /*Driver for Runge-Kutta solver*/
float ystart[],x1,x2,eps,h1,hmin;
int nvar,*nok,*nbad;
void (*derivs)(); /* ANSI: void (*derivs)(float,float *,float *); */
void (*rkqc)(); /* ANSI: void (*rkqc)(float *,float *,int,float
*,float,
float,float *,float *,float *,void (*)()); */
{
int nstp,i;
float xsav,x,hnext,hdid,h;
float *yscal,*y,*dydx,*vector();
void nrerror(),free_vector();
yscal=vector(1,nvar); /*Define arrays so that they start at 1.
This is necessary because the routines
were ported from
Fortran.*/
y=vector(1,nvar);
dydx=vector(1,nvar);
x=x1;
h=(x2 > x1) ? fabs(h1) : -fabs(h1);
*nok = (*nbad) = kount = 0;
for (i=1;i<=nvar;i++) y[i]=ystart[i];
if (kmax > 0) xsav=x-dxsav*2.0; /*Guarantee storage of first
step*/
for (nstp=1;nstp<=MAXSTP;nstp++) { /*Take MAXSTP at the most*/
(*derivs)(x,y,dydx);
for (i=1;i<=nvar;i++) /*Scale accuracy*/
yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
if (kmax > 0) {
if (fabs(x-xsav) > fabs(dxsav)) { /*Store intermediate values*/
if (kount < kmax-1) {
xp[++kount]=x;
for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
xsav=x;
}
}
}
if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x; /*Cut down step size?*/
(*rkqc)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);
if (hdid == h) ++(*nok); else ++(*nbad);
if ((x-x2)*(x2-x1) >= 0.0) { /*Are we done?*/
for (i=1;i<=nvar;i++) ystart[i]=y[i];
if (kmax) {
xp[++kount]=x;
for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
}
free_vector(dydx,1,nvar);
free_vector(y,1,nvar);
free_vector(yscal,1,nvar);
return; /*Normal exit*/
}
if (fabs(hnext) <= hmin) nrerror("Step size too small in
ODEINT");
h=hnext;
}
nrerror("Too many steps in routine ODEINT");
}
main()
{
int i,nbad,nok,j;
float eps,h1,hmin,x1,x2,*ystart;
ystart=vector(1,N);
x1=0.0; /*Initial time step*/
x2=0.002;
ystart[1]=0.0; /*Set the initial values for the variables*/
ystart[2]=0.0;
ystart[3]=1.72;
ystart[4]=100.0;
ystart[5]=10.0;
ystart[6]=0.0;
ystart[7]=0.0;
ystart[8]=0.0;
ystart[9]=0.0;
ystart[10]=0.0;
ystart[11]=3.49;
ystart[12]=4.8;
eps=1.0e-4; /*Error tolerance*/
h1=1.0;
hmin=0.0;
for(j=1;j<=ITER;j++) {
x1=x1+0.002; /*Increase time step*/
x2=x2+0.002;
odeint(ystart,N,x1,x2,eps,h1,hmin,&nok,&nbad,derivs,rkqc);
/*Call Runge-Kutta
Driver*/
for(i=1;i<=N;i++) {
printf("%f ",ystart[i]);
}
printf("\n");
}
free_vector(ystart,1,N);
}