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