sub new {
	my $arg1 = shift;
	my $class = ref($arg1) || $arg1; # can be used as class or instance method
	my $self  = {};   # ref to an empty hash
	bless $self, $class;
	$self->_initialise();
	return $self;
}

sub rk2 { my ($ynref, $dydtref, $t, $dt) = @_;
	if (ref $ynref ne 'ARRAY') {
		warn "Math::RungeKutta::rk2: 1st arg must be an arrayref\n"; return ();
	}
	if (ref $dydtref ne 'CODE') {
		warn "Math::RungeKutta::rk2: 2nd arg must be a subroutine ref\n";
		return ();
	}

	my $gamma = .75;  # Ralston's minimisation of error bounds
	my $alpha = 0.5/$gamma; my $beta = 1.0-$gamma;
	my $alphadt=$alpha*$dt; my $betadt=$beta*$dt; my $gammadt=$gamma*$dt;
	my $ny = $#$ynref;
	my @ynp1; $#ynp1 = $ny;
	my @dydtn; $#dydtn = $ny;
	my @ynpalpha; $#ynpalpha = $ny;  # Gear calls this q
	my @dydtnpalpha; $#dydtnpalpha = $ny;

	@dydtn = &{$dydtref}($t, @$ynref);
	my $i; for ($i=$[; $i<=$ny; $i++) {
		$ynpalpha[$i] = ${$ynref}[$i] + $alphadt*$dydtn[$i];
	}
	@dydtnpalpha = &{$dydtref}($t+$alphadt, @ynpalpha);
	for ($i=$[; $i<=$ny; $i++) {
		$ynp1[$i] = ${$ynref}[$i]+$betadt*$dydtn[$i]+$gammadt*$dydtnpalpha[$i];
	}
	return ($t+$dt, @ynp1);
}
my @saved_k0; my $use_saved_k0 = 0;
sub rk4 { my ($ynref, $dydtref, $t, $dt) = @_;
	# The Runge-Kutta-Merson 5-function-evaluation 4th-order method
	# in the sine-cosine example, this seems to work as a 7th-order method !
	if (ref $ynref ne 'ARRAY') {
		warn "Math::RungeKutta::rk4: 1st arg must be an arrayref\n";
		return ();
	}
	if (ref $dydtref ne 'CODE') {
		warn "Math::RungeKutta::rk4: 2nd arg must be a subroutine ref\n";
		return ();
	}
	my $ny = $#$ynref; my $i;

	# @eta0 = @#ynref;
	my @k0; $#k0=$ny;
	if ($use_saved_k0) { @k0 = @saved_k0;
	} else { @k0 = &{$dydtref}($t, @$ynref);
	}
	for ($i=$[; $i<=$ny; $i++) { $k0[$i] *= $dt; }

	my @eta1; $#eta1=$ny;
	for ($i=$[; $i<=$ny; $i++) { $eta1[$i] = ${$ynref}[$i] + $k0[$i]/3.0; }
	my @k1; $#k1=$ny;
	@k1 = &{$dydtref}($t + $dt/3.0, @eta1);
	for ($i=$[; $i<=$ny; $i++) { $k1[$i] *= $dt; }

	my @eta2; $#eta2=$ny;
	my @k2; $#k2=$ny;
	for ($i=$[; $i<=$ny; $i++) {
		$eta2[$i] = ${$ynref}[$i] + ($k0[$i]+$k1[$i])/6.0;
	}
	@k2 = &{$dydtref}($t + $dt/3.0, @eta2);
	for ($i=$[; $i<=$ny; $i++) { $k2[$i] *= $dt; }

	my @eta3; $#eta3=$ny;
	for ($i=$[; $i<=$ny; $i++) {
		$eta3[$i] = ${$ynref}[$i] + ($k0[$i]+3.0*$k2[$i])*0.125;
	}
	my @k3; $#k3=$ny;
	@k3 = &{$dydtref}($t+0.5*$dt, @eta3);
	for ($i=$[; $i<=$ny; $i++) { $k3[$i] *= $dt; }

	my @eta4; $#eta4=$ny;
	for ($i=$[; $i<=$ny; $i++) {
		$eta4[$i] = ${$ynref}[$i] + ($k0[$i]-3.0*$k2[$i]+4.0*$k3[$i])*0.5;
	}
	my @k4; $#k4=$ny;
	@k4 = &{$dydtref}($t+$dt, @eta4);
	for ($i=$[; $i<=$ny; $i++) { $k4[$i] *= $dt; }

	my @ynp1; $#ynp1 = $ny;
	for ($i=$[; $i<=$ny; $i++) {
		$ynp1[$i] = ${$ynref}[$i] + ($k0[$i]+4.0*$k3[$i]+$k4[$i])/6.0;
	}

	# Merson's method for error estimation, see Gear p85, only works
	# if F is linear, ie F = Ay + bt, so that eta4 has no 4th-order
	# errors.  So in general step-doubling is the only way to do it.
	# Estimate error terms ...
	# if ($epsilon) {
	# 	my $errmax = 0; my $diff;
	# 	for ($i=$[; $i<=$ny; $i++) {
	# 		$diff = 0.2 * abs ($ynp1[$i] - $eta4[$i]);
	# 		if ($errmax < $diff) { $errmax = $diff; }
	# 	}
	# 	# print "errmax = $errmax\n"; # not much related to the actual error
	# }

	return ($t+$dt, @ynp1);
}
my $t; my $halfdt; my @y2;
sub rk4_auto { my $ynref=shift; my $dydtref=shift; $t=shift;
	my ($dt, $arg4) = @_;
	if (ref $ynref ne 'ARRAY') {
		warn "Math::RungeKutta::rk4_auto: 1st arg must be an arrayref\n";
		return ();
	}
	if (ref $dydtref ne 'CODE') {
		warn "Math::RungeKutta::rk4_auto: 2nd arg must be a subroutine ref\n";
		return ();
	}
	if ($dt == 0) { $dt = 0.1; }
	my @errors; my $epsilon;
	if (ref $arg4 eq 'ARRAY') {
		@errors = @$arg4; undef $epsilon;
	} else {
		$epsilon = abs $arg4; undef @errors;
		if (! $epsilon) { $epsilon = .0000001; }
	}
	my $ny = $#$ynref; my $i;

	my @y1; $#y1 = $#$ynref;
	$#y2 = $#$ynref;
	my @y3; $#y3 = $#$ynref;
	$#saved_k0 = $#$ynref; @saved_k0 = &{$dydtref}($t, @$ynref);
	my $resizings = 0;
	my $highest_low_error = 0.1E-99; my $highest_low_dt = 0.0;
	my $lowest_high_error = 9.9E99;  my $lowest_high_dt = 9.9E99;
	while (1) {
		$halfdt = 0.5 * $dt; my $dummy;
		$use_saved_k0 = 1;
		($dummy, @y1) = &rk4($ynref, $dydtref, $t, $dt);
		($dummy, @y2) = &rk4($ynref, $dydtref, $t, $halfdt);
		$use_saved_k0 = 0;
		($dummy, @y3) = &rk4(\@y2, $dydtref, $t+$halfdt, $halfdt);

		my $relative_error;
		if ($epsilon) {
	 		my $errmax = 0; my $diff; my $ymax = 0;
	 		for ($i=$[; $i<=$ny; $i++) {
	 			$diff = abs ($y1[$i]-$y3[$i]);
	 			if ($errmax < $diff) { $errmax = $diff; }
	 			if ($ymax < abs ${$ynref}[$i]) { $ymax = abs ${$ynref}[$i]; }
	 		}
			$relative_error = $errmax/($epsilon*$ymax);
		} elsif (@errors) {
			$relative_error = 0.0; my $diff;
	 		for ($i=$[; $i<=$ny; $i++) {
	 			$diff = abs ($y1[$i]-$y3[$i]) / abs $errors[$i];
	 			if ($relative_error < $diff) { $relative_error = $diff; }
	 		}
                        if ($relative_error ==0) {$relative_error = abs($errors[0]);}  # TRM added this
		} else {
			die "RungeKutta::rk4_auto: \$epsilon & \@errors both undefined\n";
		}
		# Gear's "correction" assumes error is always in 5th-order terms :-(
		# $y1[$i] = (16.0*$y3{$i] - $y1[$i]) / 15.0;
		if ($relative_error < 0.60) {
			if ($dt > $highest_low_dt) {
				$highest_low_error = $relative_error; $highest_low_dt = $dt;
			}
		} elsif ($relative_error > 1.67) {
			if ($dt < $lowest_high_dt) {
				$lowest_high_error = $relative_error; $lowest_high_dt = $dt;
			}
		} else {
			last;
		}
		if ($lowest_high_dt<9.8E99 && $highest_low_dt>1.0E-99) { # interpolate
			my $denom = log ($lowest_high_error/$highest_low_error);
			if ($highest_low_dt==0.0||$highest_low_error==0.0||$denom == 0.0) {
				$dt = 0.5 * ($highest_low_dt+$lowest_high_dt);
			} else {
				$dt = $highest_low_dt * ( ($lowest_high_dt/$highest_low_dt)
				 ** ((log (1.0/$highest_low_error)) / $denom) );
			}
		} else {
			my $adjust = $relative_error**(-0.2); # hope error is 5th-order ...
			if (abs $adjust > 2.0) {
				$dt *= 2.0;  # prevent infinity if 4th-order is exact ...
			} else {
				$dt *= $adjust;
			}
		}
		$resizings++;
		if ($resizings>4 && $highest_low_dt>1.0E-99) {
			# hope a small step forward gets us out of this mess ...
			$dt = $highest_low_dt;  $halfdt = 0.5 * $dt;
			$use_saved_k0 = 1;
			($dummy, @y2) = &rk4($ynref, $dydtref, $t, $halfdt);
			$use_saved_k0 = 0;
			($dummy, @y3) = &rk4(\@y2, $dydtref, $t+$halfdt, $halfdt);
			last;
		}
	}

	return ($t+$dt, $dt, @y3);
}
sub rk4_auto_midpoint {
	return ($t+$halfdt, @y2);
}

# ---------------------- EXPORT_OK routines ----------------------

sub rk4_ralston { my ($ynref, $dydtref, $t, $dt) = @_;
	if (ref $ynref ne 'ARRAY') {
		warn "RungeKutta::rk4_ralston: 1st arg must be arrayref\n"; return ();
	}
	if (ref $dydtref ne 'CODE') {
		warn "RungeKutta::rk4_ralston: 2nd arg must be a subroutine ref\n";
		return ();
	}
	my $ny = $#$ynref; my $i;

	# Ralston's minimisation of error bounds, see Gear p36
	my $alpha1=0.4; my $alpha2 = 0.4557372542; # = .875 - .1875*(sqrt 5);

	if ($use_saved_k0) { @k0 = @saved_k0;
	} else { @k0 = &{$dydtref}($t, @$ynref);
	}
	for ($i=$[; $i<=$ny; $i++) { $k0[$i] *= $dt; }

	my @k1; $#k1=$ny;
	for ($i=$[; $i<=$ny; $i++) { $k1[$i] = ${$ynref}[$i] + 0.4*$k0[$i]; }
	@k1 = &{$dydtref}($t + $alpha1*$dt, @k1);
	for ($i=$[; $i<=$ny; $i++) { $k1[$i] *= $dt; }

	my @k2; $#k2=$ny;
	for ($i=$[; $i<=$ny; $i++) {
		$k2[$i] = ${$ynref}[$i] + 0.2969776*$k0[$i] + 0.15875966*$k1[$i];
	}
	@k2 = &{$dydtref}($t + $alpha2*$dt, @k2);
	for ($i=$[; $i<=$ny; $i++) { $k2[$i] *= $dt; }

	my @k3; $#k3=$ny;
	for ($i=$[; $i<=$ny; $i++) {
		$k3[$i] = ${$ynref}[$i] + 0.21810038*$k0[$i] - 3.0509647*$k1[$i]
		 + 3.83286432*$k2[$i];
	}
	@k3 = &{$dydtref}($t+$dt, @k3);
	for ($i=$[; $i<=$ny; $i++) { $k3[$i] *= $dt; }

	my @ynp1; $#ynp1 = $ny;
	for ($i=$[; $i<=$ny; $i++) {
		$ynp1[$i] = ${$ynref}[$i] + 0.17476028*$k0[$i]
		 - 0.55148053*$k1[$i] + 1.20553547*$k2[$i] + 0.17118478*$k3[$i];
	}
	return ($t+$dt, @ynp1);
}
sub rk4_classical { my ($ynref, $dydtref, $t, $dt) = @_;
	if (ref $ynref ne 'ARRAY') {
		warn "Math::RungeKutta::rk4_classical: 1st arg must be arrayref\n";
		return ();
	}
	if (ref $dydtref ne 'CODE') {
		warn "RungeKutta::rk4_classical: 2nd arg must be subroutine ref\n";
		return ();
	}
	my $ny = $#$ynref; my $i;

	# The Classical 4th-order Runge-Kutta Method, see Gear p35

	if ($use_saved_k0) { @k0 = @saved_k0;
	} else { @k0 = &{$dydtref}($t, @$ynref);
	}
	for ($i=$[; $i<=$ny; $i++) { $k0[$i] *= $dt; }

	my @eta1; $#eta1=$ny;
	for ($i=$[; $i<=$ny; $i++) { $eta1[$i] = ${$ynref}[$i] + 0.5*$k0[$i]; }
	my @k1; $#k1=$ny;
	@k1 = &{$dydtref}($t+0.5*$dt, @eta1);
	for ($i=$[; $i<=$ny; $i++) { $k1[$i] *= $dt; }

	my @eta2; $#eta2=$ny;
	for ($i=$[; $i<=$ny; $i++) { $eta2[$i] = ${$ynref}[$i] + 0.5*$k1[$i]; }
	my @k2; $#k2=$ny;
	@k2 = &{$dydtref}($t+0.5*$dt, @eta2);
	for ($i=$[; $i<=$ny; $i++) { $k2[$i] *= $dt; }

	my @eta3; $#eta3=$ny;
	for ($i=$[; $i<=$ny; $i++) { $eta3[$i] = ${$ynref}[$i] + $k2[$i]; }
	my @k3; $#k3=$ny;
	@k3 = &{$dydtref}($t+$dt, @eta3);
	for ($i=$[; $i<=$ny; $i++) { $k3[$i] *= $dt; }

	my @ynp1; $#ynp1 = $ny;
	for ($i=$[; $i<=$ny; $i++) {
		$ynp1[$i] = ${$ynref}[$i] +
		 ($k0[$i] + 2.0*$k1[$i] + 2.0*$k2[$i] + $k3[$i]) / 6.0;
	}
	return ($t+$dt, @ynp1);
}

# --------------------- infrastructure ----------------------

sub arr2txt { # neat printing of arrays for debug use
	my @txt; foreach (@_) { push @txt, sprintf('%g',$_); }
	return join (' ',@txt)."\n";
}
my $flag;
sub gaussn {   my $standdev = $_[$[];
	# returns normal distribution around 0.0 by the Box-Muller rules
	if (! $flag) {
		$a = sqrt(-2.0 * log(rand)); $b = 6.28318531 * rand;
		$flag = 1; return ($standdev * $a * sin($b));
	} else {
		$flag = 0; return ($standdev * $a * cos($b));
	}
}

1;