--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+expint/expm_Arnoldi.m	Mon Jun 26 20:15:54 2017 +0200
@@ -0,0 +1,54 @@
+function y = expm_Arnoldi(A,v,t,toler,m)
+%
+% y = expm_Arnoldi(A,v,t,toler,m)
+%
+% computes $y = \exp(-t A) v$
+% input: A (n x n)-matrix, v n-vector, t>0 time interval,
+%        toler>0 tolerance, m maximal Krylov dimension
+%
+% Copyright (c) 2012 by M.A. Botchev
+% Permission to copy all or part of this work is granted,
+% provided that the copies are not made or distributed
+% for resale, and that the copyright notice and this
+% notice are retained.
+%
+% THIS WORK IS PROVIDED ON AN "AS IS" BASIS.  THE AUTHOR
+% PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED,
+% REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS
+% MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE.
+%
+n = size (v,1);
+V = zeros(n  ,m+1);
+H = zeros(m+1,m);
+
+beta = norm(v);
+V(:,1) = v/beta;
+resnorm = inf;
+
+j=0;
+
+while resnorm > toler
+    j = j+1;
+    w = A*V(:,j);
+    for i=1:j
+        H(i,j) = w'*V(:,i);
+	w      = w - H(i,j)*V(:,i);
+    end
+    H(j+1,j) = norm(w);
+    e1 = zeros(j,1); e1(1) = 1;
+    ej = zeros(j,1); ej(j) = 1;
+    s  = [0.01, 1/3, 2/3, 1]*t;
+    for q=1:length(s)
+        u         = expm(-s(q)*H(1:j,1:j))*e1;
+        beta_j(q) = -H(j+1,j)* (ej'*u);
+    end
+    resnorm = norm(beta_j,'inf');
+   % fprintf('j = %d, resnorm = %.2e\n',j,resnorm);
+    if resnorm<=toler
+       break
+    elseif j==m
+       disp('warning: no convergence within m steps');
+    end
+    V(:,j+1) = w/H(j+1,j);
+end
+y = V(:,1:j)*(beta*u);