function variationaleq_ex
% Simple nonlinear autonomous example 
% dotx1 = -x1^3x2  and dotx2 = -x1x2^2, with x1(0)=x2(0)=1.
% J = Jacobian matrix, so the associated matrix variational equation is dotXI = J XI, XI(0) = eye(2).
tspan = [0,-exp(-1)+5.0e-3];
% Initial conditions for x_1, x_2, and xi = eye(2)
Y0 = [1,1,1,0,0,1];
% Integrate to reasonably tight tolerances
opts = odeset( 'reltol', 1.0e-8, 'abstol', 1.0e-8 );
% Put the solution in a structure, to compute the residual.
sol = ode45( @bothodes, tspan, Y0, opts );

% We refine the mesh ourselves so as to be sure that our residual computation reflects the actual changes in the solution as found by ode45.
nRefine = 9;
n = length( sol.x );
size( sol.x )
h = diff( sol.x );
tRefine = RefineMesh( sol.x, nRefine );
numpoints = length(tRefine);

% Now compute the residual.
[yhat,yphat] = deval(sol,tRefine);
Resid = yphat - bothodes(tRefine,yhat);
% The residual relative to the size of the rhs is also of interest.
RResid = yphat./bothodes(tRefine,yhat) - 1;

figure(1),plot(tRefine,Resid(1,:),'k-')
set(gca,'fontsize',16);
title('Residual in x_1')
figure(2),plot(tRefine,Resid(2,:),'k-')
set(gca,'fontsize',16);
title('Residual in x_2')
figure(3), semilogy(tRefine,abs(RResid), 'k-')
axis([-0.4,0,10E-12,10E-6]);
set(gca,'fontsize',16);
title('Relative Residual')

% Now look at the condition number
sigma1 = zeros(1,n);
sigma2 = zeros(1,n);
cond  = zeros(1,n);
for k=1:n,
    Xt=[sol.y(3,k),sol.y(4,k);sol.y(5,k),sol.y(6,k)];
    sigma = svd(Xt);
    sigma1(k) = sigma(1);
    sigma2(k) = sigma(2);
    cond(k) = sigma1(k)/min(sigma2(1:k));
end
figure(4), semilogy(sol.x,cond,'k')
set(gca,'fontsize',16);
title('Condition number')
end

function yp=bothodes(t,y)
    yp=zeros(size(y));
    %yp=[x1;x2;xi11;xi12;xi21;xi22];
    yp(1,:)=-y(1,:).^3.*y(2,:);
    yp(2,:)=-y(1,:).*y(2,:).^2;
    yp(3,:)=-3*y(1,:).^2.*y(2,:).*y(3,:)-y(1,:).^3.*y(5,:);
    yp(4,:)=-3*y(1,:).^2.*y(2,:).*y(4,:)-y(1,:).^3.*y(6,:);
    yp(5,:)=-y(2,:).^2.*y(3,:)-2*y(1,:).*y(2,:).*y(5,:);
    yp(6,:)=-y(2,:).^2.*y(4,:)-2*y(1,:).*y(2,:).*y(6,:);
end


