% Method of modified equations exotic example
% x'' + x = 0
%
% 2nd order Taylor series method (for x')
% exact integration of x' to get x
%
a = 0;
b = 5*pi;
n = 30;
h = (b-a)/n;
x = zeros(1,n+1);
x(1) = 1;
dx = zeros(1,n+1);
dx(1) = 0;

% Integration steps
for i=2:n+1,
    dx(i) = dx(i-1) - h*x(i-1) - h^2*dx(i-1)/2;
    x(i) = x(i-1) + h*dx(i-1) - h^2*x(i-1)/2 - h^3*dx(i-1)/6;
end
% Integration complete.

tc = linspace(a,b,n+1);
tf = linspace(a,b,5*n+1);
xref = cos( tf );
xrefc = cos( tc );

% Analytical solution of modified equation
w = 1 + h^2/12;
zet = -h^3/48;
alf = x(1);
bet = (dx(1) + zet*w*alf)/w/sqrt(1-zet^2) ; 

% fine and coarse
modxref = exp(-zet*w*tf).*(alf*cos(w*sqrt(1-zet^2)*tf) ...
                         + bet*sin(w*sqrt(1-zet^2)*tf) );
modxrefc = exp(-zet*w*tc).*(alf*cos(w*sqrt(1-zet^2)*tc) ...
                          + bet*sin(w*sqrt(1-zet^2)*tc) );
% Second plot shows the errors clearly; first plot is nicer.
figure(1),plot( tc, x, 'k+', tf, xref, 'k--', tf, modxref, 'k-.' );
set(gca,'fontsize',16);
axis([0,16,-1.5,1.5]);
set(gca,'XTick',0:2:16);
gure(2),semilogy( tc, abs(x-xrefc), 'k--', tc, abs(x-modxrefc), 'k-.')
set(gca,'fontsize',16);
axis([0,16,1E-4,1]);
set(gca,'XTick',0:2:16);
