function [t,U]=my_ode2(f,int,ua,h)
% my_ode2 - solves initial value problem for general system of 
%           ordinary differential equations u'=f(t,u). 
%           Mid-point method, also called cG(1). 
%   Syntax: 
%           [t,U]=my_ode(f,int,ua,h)
%   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
%           h   - positive number, the stepsize
%   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 ua 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 mid-point method, also called cG(1). 
%   Examples: 
%          [t,U]=my_ode2('funk3', [0, 20], [1; 0],.01);
%   See also:
%          ode23.m, ode45.m, ode15s.m 
%           

%----------------------------------------------------------------------
%   Bugs:
%
%   Date created: 2001-11-26
%          
%   Author:  Stig Larsson
%          
%
%   Change history:
%
%------------------------------+---------------------------------------+

a=int(1); b=int(2);   % the end points of the interval

tol = h^3;            % tolerance for fixed point iteration

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

while t(i)<b

  i=i+1;
  t(i)=t(i-1)+h; 
  tt=t(i-1)+h/2;      %  the mid-point 

  UU=U(:,i-1);        %  fixed point iteration 
  ee=tol+1;  
  while ee > tol 
    vv=UU;            % save the previous UU 
    UU=U(:,i-1)+feval(f, tt, ( U(:,i-1)+UU )/2)*h;
    ee=norm(vv-UU);   % the iteration error
  end
  U(:,i)=UU;

end
 
U=U';
t=t';


