--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+expint/Arnoldi.m	Mon Nov 06 11:42:15 2017 +0100
@@ -0,0 +1,37 @@
+function y = Arnoldi(A,v,t,tol,m)
+
+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 > tol
+    j = j+1;
+    z = A*V(:,j);
+    for i=1:j
+        H(i,j) = z'*V(:,i);
+        z  = z - H(i,j)*V(:,i);
+    end
+    H(j+1,j) = norm(z);
+    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) = z/H(j+1,j);
+end
+y = V(:,1:j)*(beta*u);