Project 6

Solve y'(t)=z(t), z'(t)=1000(1-y(t)^2)z(t)-y(t) with y(0)=2, z(0)=0 for t between 0 and 3000, using EPSODE. This problem has internal boundary layers, and the initial condition is nearly on the limit cycle, which has a perior of about 1615.5. Examine how the numerical solution depends on the choice of error tolerance.
This program uses epsode and d1mach to solve the ODE, using a stepsize h = .1 and a user-input value of epsilon. It creates three output files, named outputy, outputz, and outputyz, that contain y vs t, z vs t, and z vs y, respectively.
To compile the code, download all three of the above programs and run the commands

gfortran -g -c -o d1mach.o d1mach.f
gfortran -g -c -o epsode.o epsode.f
g++ -c hw6.C
gfortran -o hw6 hw6.o epsode.o d1mach.o -lstdc++

the code is then executed by the command

hw6

The variable "index" controls how epsode solves the ODE. In the current implementation, "index" is set to 22, which has epsode use the pederv function to evaluate the Jacobian exactly in solving the ODE. With this implementation, the solution converges at eps = 10^-8. Graphs of the solution at various values of epsilon are below.

Z vs Y, eps = 10^-5
Z vs Y, eps = 10^-6
Z vs Y, eps = 10^-7
Z vs Y, eps = 10^-8
Z vs Y, eps = 10^-9

Z vs t, eps = 10^-8
Z vs t, eps = 10^-9
Y vs t, eps = 10^-8
Y vs t, eps = 10^-9