% [tout, yout] = chebode(f, y0, t, n)
%
% Solve a scalar ODE y' = f(t,y) by n steps of a Runge-Kutta-Chebyshev
% method, with t = [t0, tfinal] and y(t0) = y0. Assume the function f
% has the signature
%
% [f_out,dfdy_out] = f(t,y)
%
% where f_out is a vector of function values and dfdy_out is a vector
% of derivatives.

function [tout, yout] = chebode(f, y0, t, n)

h = (t(2)-t(1))/n;
[D,x] = cheb(8);
tout = linspace(t(1), t(2), n+1)';
yout = 0*tout;
yout(1) = y0;

Ds = 2*D/h;
v = Ds(1:end-1,end);
Ds = Ds(1:end-1,1:end-1);
xs = h*(1+x)/2;
xs = xs(1:end-1);

for k = 1:n
tc = tout(k) + xs;

% Approximate solution by a constant over the step
yc = 0*tc + yout(k);

% Get initial residual information
[fc,fpc] = f(tc,yc);
r = Ds*yc-fc+v*yout(k);

for j = 1:20

% Compute Newton direction
dyc = ( Ds-diag(fpc) )\r;

% Line search
yc_old = yc;
normr_old = norm(r);
s = 1;
for i = 1:20
yc = yc_old-s*dyc;
[fc,fpc] = f(tc,yc);
r = Ds*yc-fc+v*yout(k);
if norm(r) < normr_old, break; end
s = s/2;
end

% Check for termination
if norm(dyc) < 1e-8*norm(yc), break; end

end

yout(k+1) = yc(1);
end

% CHEB compute D = differentiation matrix, x = Chebyshev grid
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D')); % diagonal entries