function [t,x] = irk3(Ffile,DFfile,t0,tf,x0,h)
% function [t,x] = irk3(Ffile,DFfile,t0,tf,x0,h)                 4/23/92
% Integrates ODE's backwards with a third order implicit Runge-Kutta method
% (implicit equations are solved via linear approximation)
%
% For the system of ODE's: dx/dt = F(x)
% IRK3.M iterates:        x[k+1] = x[k] -h [g G1 + G2 J1^(-1) J2]^(-1) F(x[k])
% G1, G2, J1, J2 are functions of DF(x[k])
% starting from x[0] = x0  (x is a column vector)
% See Mulligan, Casey B.  "Progressive Paths" equation (6-4)
%
% INPUTS:
% Ffile  = filename of F(x), a string
% DFfile = filename of DF(x), a string
% Ffile and DFfile should be function files with one (vector) argument, and
% Ffile should have column vector output of the same dimension as x
% DFfile should output a square matrix, a row the partials of one component
% Integration is from t0 to tf, at step size h
% OUTPUTS:
% t is a row vector of times [t0 (t0+h) (t0+2h) ... tf]
% x is a matrix of points on the trajectory [x(t0) x(t0+h) x(t0+2h) ...x(tf)]

% Implicit Two Stage Runge-Kutta parameters
beta11 = 1/2;  beta12 = 2/3;  beta21 = -1/8;  beta22 = 1/2;  gamma = 3/19;

xk = x0; tk = t0; x = xk; t = t0; n = length(x0);
while tk <= tf
     eval(['Fofx = ' Ffile '(xk);']);
     eval(['DFofx = ' DFfile '(xk);']);
     G1 = eye(n)/(eye(n) + beta11*h*DFofx);
     G2 = (1-gamma)*eye(n) - beta12*gamma*G1*h*DFofx;
     J1 = eye(n) + beta22*h*DFofx - beta12*beta21*h*DFofx*G1*h*DFofx;
     J2 = eye(n) - beta21*h*DFofx*G1;
     xk = xk - h*( gamma*G1 + G2*(J1\J2) )*Fofx;
     tk = tk + h;
     t = [t tk];  x = [x xk];
end
