function twostep_test1

% CCK 3018-Mar-10 minor debugging to run with Matlab.

% Initial conditions
x0 = [1,0];  % State is expected to be a row vector.
t0 = 0;

% Set up for the integration with coarse steps
tfinal = 2*pi;
Nsteps = 10;


%%%%%%%%%%%%%%%%%% This part is worth looking at closely. %%%%%%%%%%%%%%%%

% Setup for twostep:
stepper = @euler;
order = 1;
solver = @(ODE,x0,t0,h) twostep(ODE,x0,t0,h,stepper,order);

% Integrate
[xa,ta] = ode_fixed(solver, @Harmonic, x0, t0, tfinal, Nsteps);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Setup for twostep:
stepper = @rk2;
order = 2;
solver = @(ODE,x0,t0,h) twostep(ODE,x0,t0,h,stepper,order);

% Integrate
[xb,tb] = ode_fixed(solver, @Harmonic, x0, t0, tfinal, Nsteps);


% Plots
tc = (0:1000)*2*pi/1000; % time axis for smooth plot of analytic solution
figure(1)
plot(xa(:,1), xa(:,2), 'r-x',...
     xb(:,1), xb(:,2), 'g-x',...
     cos(tc), sin(tc), 'b-' )
title('twostep\_test1.m: test of step doubling','fontsize',18)
legend('10 steps doubled euler (=rk2)','10 steps doubled rk2 (=rk3)','Analytic soln')
legend('location','NorthEast')
xlabel('x(1)','fontsize',18)
ylabel('x(2)','fontsize',18)
axis([-1.2,2,-1.2,1.2])
axis equal
% hold on
% for n = 1:numel(ta)
%    plot([0, xa(n, 1)], [0, xa(n,2)], 'r-')
% end
% hold off

figure(2)
plot (ta, xa(:,2), 'r-x', tb, xb(:,2), 'g-x', tc, sin(tc), 'b')
title('twostep\_test1.m: test of step doubling','fontsize',18)
legend('10 steps doubled euler (=rk2)','10 steps doubled rk2 (=rk3)','Analytic soln')
xlabel('t','fontsize',18)
ylabel('x(2)','fontsize',18)
legend('location','southwest')