
function [t, U] = my_ode(f, int, ua, h)

% my_ode  - solves initial value problem for general system of
%           ordinary differential equations u'=f(t,u).
%           Mid-point method, also called cG(1) (cG:continuous Galerkin method).
%   Syntax:
%           [t,U]=my_ode(f,int,ua,h)
%   Arguments:
%           f   - a function handle that evaluates the right side of the
%                 differential equations.
%           int - 1x2 matrix specifying a time interval int=[a,b]
%           ua  - dx1 matrix specifying an initial value
%           h   - positive number, the stepsize
%   Returns:
%           t - row vector containg the time points with t(1)=a
%           U - dxn 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 functionhandle
%           @f must call the function f(t,u) with syntax uprime=f(t,u).
%           The program returns the row vector t of time points with
%           t(1)=a and the dxn matrix U with column number i containing
%           the solution vector at time t(i).
%           The program uses the mid-point method, also called cG(1).
%   Examples:
%          [t,U] = my_ode(@funk, [0, 10], [0; 1],.01);
%
%

%----------------------------------------------------------------------
%
%
%   Date created: 2007-01-26
%
%   Author:  Nils Svanstedt
%   Revised: Jan Södersten      2008-01-25
%            Niklas Eriksen     2008-12-19
%
%------------------------------+---------------------------------------+

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

t = int(1):h:int(2);  % time vector   

U(:, 1) = ua;           % the first column vector of solution U
V       = ua;

for i = 2:length(t)
  tt = t(i) - h/2;      %  the mid-point
  ee = tol + 1;

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

  U(:,i) = V;
end

