This solution makes use of three MATLAB features, the integral function that does integrals adaptively and numerically, anonymous functions, which allow you to pass user-defined functions into other functions, such as integral, and the arrayfun function, which evaluates the anonymous function fk at each point of the array x.

This is Patrick Grice's code, with some editing and live-script annotation by Prof. Goodman.

By using a Live Script, I have been able to explain my code a little bit better and include the outputs in a nice little report. I have chosen to export it to HTML which looks better. PDF export is also available, but it doesn't always break the pages in nice places.

I chose five values of λ: the critical value for convergence , two above it, and two below it.

lambdas = [.2 .5 2/pi .75 2];

nLambda=length(lambdas);

f0s = {

@(x) 0, ...

@(x) 1, ...

@(x) x, ...

@(x) x.^2, ...

@(x) 10*exp(-4*x.^2) ...

};

f0string = {'$f_0(x)=0$',...

'$f_0(x)=1$',...

'$f_0(x)=x$',...

'$f_0(x)=x^2$',...

'$f_0(x)=10 e^{-4 x^2}$'};

nf0=length(f0s);

nIterates=4;

x = linspace(-1, 1, 101);

tiledlayout(length(lambdas), nf0, 'TileSpacing', 'tight', 'Padding', 'tight', ...

'TileIndexing', 'columnmajor');

for i = 1:length(f0s) % Loop over the initial functions

f0=f0s{i}; % Loop over the lambda values

for j = 1:nLambda

lambda=lambdas(j);

nexttile

fk=f0;

plot(x, arrayfun(fk, x));

hold on;

% Iterate the operator T

for k = 1:nIterates

fk = T(fk, lambda);

plot(x, arrayfun(fk, x));

end

if j == 1; title(f0string{i});end

if i == 1; ylabel(sprintf('$\\lambda=%0.2f$',lambdas(j)));end

end

end

It's very hard to see the divergence for large λ and or .

Define the operator T

function g = T(f, lambda)

g = @(x) sin(2*pi*x) + lambda * integral(@(y) f(y) ./ (1 + (x-y).^2), -1, 1);

end