program main
	use rk
	use pc
	implicit none
	integer, parameter :: N = 10000
	real, parameter :: acc = 0.001, eps = 1e-5
	real, parameter :: a = 0., b = 10.
	real :: xlist(N), ylist(N,2), h
	real :: solution, t1, t2
	integer :: i, nrk = 1
	
	xlist = 0
	xlist(1) = a
	ylist(1,1:2) = (/ 1., 1. /)
	h = 0.01
	
	! Call the driver:
	call cpu_time(t1)
	call rkdriver(xlist, ylist, b, acc, eps, h, nrk)
	call cpu_time(t2)

	print *, "*******************************"
	print *, "Runge-Kutta Method:"
	print *, ""
	print *, "Number of points:", nrk
	print *, "Elapsed time:", (t2-t1)*1000, " msec."

	! Write solutions to file:
	open(51, file='points_rk.dat', status='replace', action='write')
	do i = 1,nrk
		write(51, '(3f20.6)') xlist(i), ylist(i,1), ylist(i,2)
	enddo
	close(51)

	!=============================================================

	! Reset parameters:
	xlist(:) = 0; ylist(:,:) = 0; nrk = 1;
	xlist(1) = a
	ylist(1,1:2) = (/ 1., 1. /)
	h = 0.01

	! Call the driver with same parameters:
	call cpu_time(t1)
	call pcdriver(xlist, ylist, b, acc, eps, h, nrk)
	call cpu_time(t2)

	print *, "*******************************"
	print *, "Predictor-Corrector Method:"
	print *, ""
	print *, "Number of points:", nrk
	print *, "Elapsed time:", (t2-t1)*1000, " msec."

	! Write solutions to file:
	open(51, file='points_pc.dat', status='replace', action='write')
	do i = 1,nrk
		write(51, '(3f20.6)') xlist(i), ylist(i,1), ylist(i,2)
	enddo
	close(51)
	
	print *, "*******************************"
	print *, "Run 'make plot' to produce plots."
end program main

! Differential-equation system to solve.
! Should give a sine and cosine function.
function f(x,y)
	real, intent(in) :: x, y(:)
	real :: f(size(y))
	
	f(1) = y(2)
	f(2) = -y(1)
end function