function soln = vie_trap(N_h,T,fcn_g,fcn_k)

% This solves the integral equation
%                   t
%    Y(t) = g(t) + Int k(t,s,Y(s))ds
%                   0 
% INPUT: 
% N_h       The number of subdivisions of [0,T].
% T         [0,T] is the interval for the solution function.
% fcn_g     The handle of the driver function g(t).
% fcn_k     The handle of the kernel function k(t,s,u).
% OUTPUT:   soln is a structure with the following components.
% soln.t    The grid points at which the solution Y(t) is approximated.
% soln.y    The approximation of Y(t) at the grid points.

% The implicit trapezoidal equation is solved by simple fixed point
% iteration at each grid point in t.  For simplicity, the program uses a 
% crude means of controlling the iteration. The iteration is executed a 
% fixed number of times, controlled by the variable 'loop'.
loop = 10;  % This is much more than is usually needed.

h = T/N_h; t = linspace(0,T,N_h+1);
g_vec = fcn_g(t); 
y_vec = zeros(size(t)); y_vec(1) = g_vec(1);

for n=1:N_h
    y_vec(n+1) = y_vec(n);  % The initial estimate for the iteration.
    k_vec = fcn_k(t(n+1),t(1:n+1),y_vec(1:n+1));
    for j=1:loop  
        y_vec(n+1) = g_vec(n+1) + h*(sum(k_vec(2:n)) + ...
            (k_vec(1) + k_vec(n+1))/2);
        k_vec(n+1) = fcn_k(t(n+1),t(n+1),y_vec(n+1));
    end
end

soln.t = t;
soln.y = y_vec;

end % vie_trap