% Solving the Van der Pol equation using Backward Euler

lambda = 0.01;

dt = 0.01;
% dt = 0.001; % Try seeing what happens if you decrease the time step (this
              % may take a bit of time to run)
t = 0:dt:20;

x = NaN(1,length(t));
v = NaN(1,length(t));

% Set the initial condition
x(1) = 1.5;
v(1) = sqrt(4-1.5^2);

% Iterate using Backward Euler
for n = 1:(length(t)-1)
   y1 = x(n); 
   y2 = v(n);
      
   % Find next iterate using Newton's Method
   z = [y1; y2];
   for m = 1:1000 % Rough estimate of steps until 'converging'
       z1 = z(1);
       z2 = z(2);
       
       F = [z2 ; lambda*(1-z1^2)*z2 - z1]; % RHS of Van der Pol system
       JF = [0, 1; -2*lambda*z1*z2-1 , lambda*(1-z1^2)]; % Jacobian of system
       
       G = [y1; y2] + dt*F - z;
       JG = dt*JF - eye(2);
       z = z - JG\G; % Using mldivide
   end
   
   x(n+1) = z(1);
   v(n+1) = z(2);
end

figure, hold on
set(gca,'FontSize',12)
plot(t,x,'ko')
plot(t,v,'k-')
xlabel('t')

% Compare results to a higher-order method:
% There is a convenient discussion of numerical solutions for the Van der
% Pol equation at https://www.mathworks.com/help/matlab/examples/differential-equations.html
% I've copied some code from there, here:
tspan = [0, 20];
y0 = [sqrt(4-1.5^2); 1.5];
Mu = 0.01;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);

plot(t,y(:,1),'b-')