function [t,U]=ode_theta(f,int,u0,h,theta)
% ode_theta - solves initial value problem for general system of 
%          ordinary differential equations u'=f(t,u) 
%   Syntax: 
%           [t,U]=ode_theta(f,int,u0,h,theta)
%   Arguments: 
%           f   - string containing the name of a function file  
%           int - 1x2 matrix specifying a time interval int=[a,b]
%           u0  - dx1 matrix specifying an initial value
%           h   - positive number, the stepsize
%           theta - a number in the interval [0,1], the theta parameter
%   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)=u0.  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 theta method. 
%           theta=0    is Euler forward 
%           theta=0.5  is the midpoint method, cG(1)
%           theta=0    is Euler backward, dG(0)
%   Examples: 
%          [t,U]=ode_theta('funk3', [0, 2], [0;1], .01, .5);
%            is the midpoint method, alias cG(1)
%   See also:
%          crash.m, ode23.m, ode45.m, ode15s.m 
%           

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

tol=1e-6;   %tolerance for newton iteration
a=int(1); b=int(2);
N=ceil((b-a)/h); 

t(1)=a;
U(:,1)=u0;
i=1;

while t(i)<b
   i=i+1;
   t(i)=a+(i-1)*h;
   U(:,i)=odenewton(f,h,t(i-1),U(:,i-1),tol,theta);
end
t=t';
U=U'; 

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

function v = odenewton(f,h,t,u,tol,theta)
v = u;
w = tol + 1;
I = eye(size(u*u'));            % unit matrix dxd
tt = t + theta*h;

while norm(w)>tol
   vv = theta*v+(1-theta)*u;
   A = theta*jacobi(f,tt,vv);    % evaluate the Jacobian 
   b = -(v-h*feval(f,tt,vv)-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-6;                      % 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

 
