close all; clear; clc;

%%
mu = 4;
dydt = @(t,y)[y(2);mu*(1-y(1)^2)*y(2)-y(1)];
%dydt = @(t,y)[0 1;-1 0]*y(:) + [0; 1]*mu*(1-y(1)^2)*y(2);
k = 0.1; t = (0:round(15./k)).*k; % Time axis
y0 = [1 1]; % Initial condition

% Forwards Euler method - explicit method
y1 = zeros(numel(t),2); y1(1,:) = y0;
for ti = 2:numel(t)
    dy = k.*dydt(t(ti),y1(ti-1,:))'; % Change in y is the current gradient
    y1(ti,:) = y1(ti-1,:)+dy; % Go forward in time
end
disp('done fe');

% Backwards Euler method - implicit method
y2 = zeros(numel(t),2); y2(1,:) = y0;
for ti = 2:numel(t)
    % Fixed point iteration to find change in y
    dy = 0; j = 0; % Initial guess for dy and iteration count
    while true
        dyn = k*dydt(t(ti),y2(ti-1,:)+dy)'; % New guess
        if all(abs(dyn-dy)<1e-6) || j>100; break, end % Stopping condition
        dy = dyn; j = j+1;
    end
    yf = y2(ti-1,:)+k.*dydt(t(ti),y2(ti-1,:))'; % Go forward in time
    %dy = k*dydt(t(ti),yf)';
    y2(ti,:) = y2(ti-1,:) + dy; % Go forward in time
end
disp('done be');

% Variable timestep
y3 = y0; t3 = 0;
while t3(end)<t(end)
    % Drop the timestep until dy is small enough
    dt = k;
    while true
        dy = dt.*dydt(t(end),y3(end,:))'; % Change in y
        if all(dy<1e-2) || dt<=1e-2, break, end
        dt = dt/2; % Drop the timestep
    end
    y3(end+1,:) = y3(end,:)+dy; % Go forward in time
    t3(end+1) = t3(end) + dt; % Record the time too
end
disp('done vt');

% Crank-Nicolson 
y4 = zeros(numel(t),2);
y4(1,:) = y0;
for ti = 2:numel(t)
    dy_fe = k.*dydt(t(ti),y4(ti-1,:))'; % Change in y is the current gradient
    
    yf = y4(ti-1,:)+dy_fe; % Go forward in time
    dy_be = k*dydt(t(ti),yf)';
    
    dy = (dy_fe+dy_be)/2;
    y4(ti,:) = y4(ti-1,:)+dy; % Go forward in time
end
disp('done cn');


% 4'th order Runge?Kutta
[te,ye] = ode45(dydt,[0 t(end)],y0);
disp('done 4rk');

%% Plot
cm = lines;
figure(1); clf;
subplot(1,2,1); hold on;
plot(te,ye(:,1),'k-','linewidth',1);
plot(t,y1(:,1),'.-','color',cm(1,:),'MarkerSize',10);
plot(t,y2(:,1),'.-','color',cm(2,:),'MarkerSize',10);
plot(t3,y3(:,1),'.-','color',cm(3,:),'MarkerSize',10);
plot(t,y4(:,1),'.-','color',cm(5,:),'MarkerSize',10);
ylim([-1 1]*3);
xlim([min(t) max(t)]);
legend('4RK','FE','BE','VT','CN','location','northwest')
grid on; box on;
axis square;
xlabel('time'); ylabel('y');

subplot(1,2,2); hold on;
plot(y1(:,1),y1(:,2),'.-','color',cm(1,:),'MarkerSize',10);
plot(y2(1:90,1),y2(1:90,2),'.-','color',cm(2,:),'MarkerSize',10);
x = interp1(t3,y3(:,1),t);
y = interp1(t3,y3(:,2),t);
plot(x(1:100),y(1:100),'.-','color',cm(3,:),'MarkerSize',10);
plot(y4(1:120,1),y4(1:120,2),'.-','color',cm(5,:),'MarkerSize',10);
plot(ye(1:220,1),ye(1:220,2),'k-','linewidth',1);
ylim([-1 1]*8);
xlim([-1 1]*3);
grid on; box on;
axis square;
xlabel('y'); ylabel('y''');