function [T,W] = SolveODE(n, u0, times, h, a, d, method);

global exponent;

% One option would be to solve this exactly since we know it is linear dw/dt=A*x+g(t)
% [V,D] = eig(A);
% expA = V*diag(exp(t*diag(D)))/V
% This would not work for upwinding, for example, but one could try expm

Aw_plus_g = @(t,w) AdvDiff(t, w, h, a, d, method);

% Read ode solver documentation!

% Specifying tspan with more than two elements
% has little effect on the efficiency of computation
% but for large systems, affects memory management.

% For stiff problems, specifying the Jacobian matrix using odeset is particularly important.
% Use the Jacobian, JPattern, or Vectorized options of odeset 
% to specify information about the Jacobian. 
% If you do not supply the Jacobian then the solver estimates it numerically 
% using finite differences.

% Jacobian sparsity pattern, specified as a sparse matrix. 
% The sparse matrix contains 1s where there might be nonzero entries in the Jacobian.

% ode15s for stiff problems (will it use sparse Jacobian???)
% ode45 for slow but accurate explicit integration
tol = 1e-9; % Not too tight of a tolerance
if(0) % Explicit - way too slow!

   options = odeset('RelTol',tol);
   tic
   [T,W] = ode45(Aw_plus_g, times, u0, options);
   toc

else % Implicit -- much faster

   % Construct the sparsity pattern of the Jacobian:
   if(method==0) % Upwind
      B = ones(n,2);
      A = spdiags(B, -1:0, n, n);
   elseif(method==1) % Centered
      B = ones(n,3);
      A = spdiags(B, -1:1, n, n);
   else % Third-order upwind biased
      B = ones(n,4);
      A = spdiags(B, -2:1, n, n);
   end

   %options = odeset('RelTol',tol); % Check how much sparsity pattern helps (factor of 3)
   options = odeset('RelTol',tol,'JPattern',A);
   tic
   [T,W] = ode15s(Aw_plus_g, times, u0, options);
   toc

end   

end
