function y = expm_Arnoldi(A,v,t,toler,m)

global iter 

% 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;
iter = 0;


while resnorm > toler
    iter = iter +1;
    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);