%  This is Xmrke.m
%  Example of the use of Rke and Yvalue in MATLAB.
%  REQUIRES Rke.m, Yvalue.m, and X2rkef.m.

%  Solves van der Pol's equation on the interval [0,20] and plots the 
%  answers obtained at each step as well as several values obtained by 
%  interpolation in each step.  In a phase plane plot, the extra values 
%  are needed for a smooth graph.  

%  Set initial values and tolerances.
x = 0;
lastx = 20;
y = [1; 1];
tol = 1e-5;
threshold = [1e-5; 1e-5];

%  Initialize output arrays to initial values.    
X = 0;
Y = y';

h = 0.1;
nasteps = 0;
k0 = zeros(2,1);
while x <= lastx
  [x,y,k0,h,flag,step,ycoeff,nasteps] = ...
    Rke('X2rkef',x,y,k0,h,tol,threshold,nasteps);
  if flag ~= 0
    error(['flag = ',int2str(flag),' at x = ',num2str(x),'.']);
  end

  %  Compute several approximations in [x-step,x] to get a smooth graph.
  %  Do not produce approximations past lastx:
  for i = 3:-1:0
    xi = min(x,lastx) - i*(step/4);
    z = Yvalue(xi,x,step,ycoeff);
    X = [X; xi];
    Y = [Y; z'];
  end
end

plot(X,Y)   
legend('y(x)','y''(x)')
title('van der Pol equation')

figure
plot(Y(:,1),Y(:,2))
title('phase plane plot')
xlabel('y(x)')
ylabel('y''(x)')
   
