f = @(x,y) [y(2,:); 9*y(2,:)+10*y(1,:)];
x = 0:10;
solinit = bvpinit( x, [1,0] );
sol = bvp4c( f, @(ya,yb) [ya(1)-1; yb(1)-1], solinit );
xi = RefineMesh( sol.x, 20 );
[y,dy] = deval( sol, xi );
res = dy - f(xi,y);
size( sol.x )
% Yields [1 49]: mesh with 48 subintervals; bvp5c is better
figure(1), plot( sol.x, sol.y(1,:), 'ko', xi, y(1,:), 'k-' )
set(gca,'fontsize',16)
xlabel( 't','fontsize',16 ), ylabel( 'y','fontsize',16 )
figure(2), semilogy( sol.x(2:end), diff(sol.x), 'ko')
set(gca,'fontsize',16)
xlabel( 't','fontsize',16 ), ylabel( 'stepsizes','fontsize',16 )
figure(3), semilogy( xi, abs(res), 'k.' )
set(gca,'fontsize',16)
xlabel( 't','fontsize',16 ), ylabel( 'residual','fontsize',16 )
axis( [ 0, 10, 1.0e-8, 1] )
