T      = 10;
tol    = 10^-6;
Nmax   = 1000;
Nfixed = 69;

function M = A (t)
  M = [0, 5/(1+t); -10/(t+1), -t^2];
end

tt = linspace (0, T, 1000);
lambda = [];

for j = 1:length(tt);
  lambda = [lambda eig (A(tt(j)))];
end

title('Eigenvalues of A(t)');
xlabel('t');
ylabel('Real and imaginary part of each eigenvalue');

plot (tt, real(lambda(1,:)), ';Re(lambda1);',
      tt, imag(lambda(1,:)), ';Im(lambda1);',
      tt, real(lambda(2,:)), ';Re(lambda2);',
      tt, imag(lambda(2,:)), ';Im(lambda2);')

gset term postscript eps monochrome;
gset output "eigenvalues.eps"
replot
gset term x11

function y = f(t,x)
  y = A(t) * x;
end

[y1,t1] = ode_bdf4_var ('f', [0,T], [0; 2], tol, Nmax);
[y2,t2] = ode_rk54 ('f', [0,T], [0; 2], tol, Nmax);

y3 = ode_rk4 ('f', [0,T], [0 ;2], Nfixed);
y31 = y3(1,:);

for j = 1:(Nfixed+1)
  if ((j>length(y31)||(abs(y31(j))>10)))
     y31(j) = NaN;
  end
end

y4 = ode_bdf4_broyden ('f', [0,T], [0 ;2], Nfixed);

pause;

title('Comparison of four different ODE Solvers');
xlabel('t');
ylabel('y1(t)');

plot (t1, y1(1,:), '-@;Scheme A;', ...
      t2, y2(1,:)+1*ones(1,length(t2)), '-@;Scheme B;',
      linspace(0, T, Nfixed+1), y31+2*ones(1,Nfixed+1), '-@;Scheme C;',
      linspace(0, T, Nfixed+1), y4(1,:)+3*ones(1,Nfixed+1), '-@;Scheme D;');

gset term postscript eps monochrome;
gset output "solution.eps"
replot
gset term x11
