function [t,U]=ode_dg0(f,int,ua,tol)
% ode_dg0 - solves initial value problem for general system of 
%          ordinary differential equations u'=f(t,u) 
%   Syntax: 
%           [t,U]=ode_dg0(f,int,ua,tol)
%   Arguments: 
%           f   - string containing the name of a function file  
%           int - 1x2 matrix specifying a time interval int=[a,b]
%           ua  - dx1 matrix specifying an initial value
%           tol   - positive number, the local tolerance
%   Returns: 
%           t - nx1 matrix containg the time points with t(1)=a
%           U - nxd matrix containing the approximate solution
%   Description:
%           The program computes an approximate solution of the intial 
%           value problem  u'=f(t,u), a<t<b; u(a)=ua.  Here u, f, 
%           and u0 are column vectors of dimension d.  The file f.m 
%           must contain the function f(t,u) with syntax uprime=f(t,u). 
%           The program returns the nx1 matrix t of time points with 
%           t(1)=a and the nxd matrix U with row number i containing 
%           the transposed solution vector at time t(i).  
%           The program uses the Euler backward method, alias dG(0),
%           with adaptivity. 
%   Examples: 
%          [t,U]=ode_dg0('funk3', [0, 2], [0; 1],.01);
%   See also:
%          crash.m, ode23.m, ode45.m, ode15s.m 
%           

%----------------------------------------------------------------------
%   Bugs:
%
%   Date created: 
%          
%   Author:
%          
%
%   Change history:
%
%------------------------------+---------------------------------------+

c=2; hmax=1;          % parameters for step selection
a=int(1); b=int(2);   % the end points of the interval

ntol=tol;
h=tol;
t(1)=a;
U(:,1)=ua;            % the first column vector of U
i=1;

while t(i)<b
  i=i+1;
  t(i)=t(i-1)+h;
  U(:,i)=odenewton(f,h,t(i-1),U(:,i-1),ntol);
  residual=norm(U(:,i)-U(:,i-1),inf)/h;         % compute residual
  error=h*residual;                             % error bound
  hh=tol/(residual+eps);
  hh=min([hh c*h hmax]);
  h=max([hh h/c]);                              % compute new step
  if error>tol                   % reject if error bound too large
      i=i-1;
  end
end

t=t';
U=U'; 

%--------------------------------------------------------------------
% subfunctions 

function v = odenewton(f,h,t,u,tol)
v = u;
w = tol + 1;
I = eye(size(u*u'));           % unit matrix dxd
s = t+h;
while norm(w)>tol
   A = jacobi(f,s,v);          % evaluate the Jacobian A=Df(t,v)
   b = -(v-h*feval(f,s,v)-u);  % evaluate the residual 
   w = (I-h*A)\b;              % solve the linearized equation
   v = v + w;                  % update 
end

function A=jacobi(f,t,x)
h = 1e-8;                      % stepsize 
I = eye(size(x*x'));           % unit matrix dxd

for i = 1:size(x,1)
   % I(:,i) is the i:th basis vector
   % compute i:th column of A 
   A(:,i) = ( feval(f,t,x+h*I(:,i)) - feval(f,t,x-h*I(:,i)) )/(2*h);
end

