Although this question is more than a decade old, there still seems to be some interest in this. I was suprised to see that all the answers use the solution y = exp(x.^2 -1) which is what we are supposed to find by solving the ODE dy/dx = 2xy, y(1) = 1 numerically.
clear
clc
% Solve: dy/dx = 2xy, y(x), y(1) = 1
N = 20;
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N); %dom=[-1,1]
% Domain is important as we will be solving for y only inside this domain.
% The ODE problem can now be written as the following linear problem with
% one constraint. This constraint arises due to the boundary condition.
% D*y_est == 2*x.*y_est
% y(x=1) = 1.
y_est = nan(N+1,1);
y_est(end) = 1; % boundary condition
% Note that we have N unknowns (values of y at x(1:N)) but the system above has N+1 equations.
% To solve this we utilize the boundary condition (constraint) to reduce
% the system to the form q = Ap with size(A) = [N,N] (N equations).
% Solving this using mldivide (\ operator) gives us
y_est(1:N) = (D(1:N,1:N)-2*diag(x(1:N)))\(-D(1:N,N+1)*1);
% Remark that we didn't use y = exp(x.^2 -1). In the figure below we can
% see that the solution to the ODE we calculated matches the exact solution
% y = exp(x.^2 -1) well at the x-grid points.
% Plotting
plot(x,y_est,'s',x,exp(x.^2-1))
xlabel("x")
ylabel("y, y_{est}")
legend('Chebyshev','Exact')
%%
% Checking convergence. Here we can see how the numerically computed
% solution improves as we increase number of x-grid points N.
figure
for N = 2:30
% compute the chebyshev differentiation matrix and x-grid
[D,x] = cheb(N); %dom=[-1,1]
% Domain is important as we will be solving for y only inside this domain.
% The ODE problem can now be written as the following linear problem with
% one constraint. This constraint arises due to the boundary condition
% y(x=1) = 1.
% D*y_est == 2*x.*y_est
y_est = nan(N+1,1);
y_est(end) = 1; % boundary condition
% Note that we have N unknowns (values of y at x(1:N)) but the system above has N+1 equations.
% To solve this we utilize the boundary condition (constraint) to reduce
% the system to the form q = Ap with size(A) = [N,N].
% Solving this using mldivide (\ operator) gives us
y_est(1:N) = (D(1:N,1:N)-2*diag(x(1:N)))\(-D(1:N,N+1)*1);
semilogy(N,norm(y_est-exp(x.^2-1)),'bo')
xlabel('N')
ylabel("||y-y_{est}||")
hold on
drawnow
end
function [D,x] = cheb(N)
% check the base case
if N == 0; D = 0; x = 1; return; end
% create the Chebyshev grid
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)));
D = D - diag(sum(D'));
end
