% Difference methods: euler's, heun, trapezoid, integration.
% Representation of 1st order differential equations (DEs) (phase portrait)
% Solve for specific intital value using differnece methods and ode45().
% Solve system of 1st order DEs.
% Solve higher order DEs as system of 1st order DEs.
% Boundary value problems using finite-differences and linear system of equations (sparse, diagonal matrix)
%% 1.

f = @(t, u) cos(3*t) - sin(5*t).*u;
a = 0;
b = 5;
u0 = 2;

N = 20;
t = linspace(a,b,N);
u = linspace(0,3,N);

[T,U] = meshgrid(t,u);
DT = ones(size(T));
DU = f(T, U);

quiver(T,U,DT,DU)
xlim([0,5]); ylim([0,3])

%
%% Uppgift 2.

close all

f = @(t, u) cos(3*t) - sin(5*t).*u;

a = 0;
b = 5;
u0 = 2;

N = 20;
t = linspace(a,b,N);
u = linspace(0,3,N);

[T,U] = meshgrid(t,u);
DT = ones(size(T));
DU = f(T, U);

quiver(T,U,DT,DU)
xlim([0,5]); ylim([0,3])
hold on


N = 200;
t = linspace(a,b,N);
h = (b-a)/N;

euler_forward_step = @(t, u, h) u + h*f(t,u);

u = zeros(length(N));
u(1) = u0;
for k = 2:N
    u(k) = euler_forward_step(t(k), u(k-1), h);
end
euler_u = u;
plot(t,u)

%
%% 3.

close all

f = @(t, u) cos(3*t) - sin(5*t).*u;

a = 0;
b = 5;
u0 = 2;

N = 20;
t = linspace(a,b,N);
u = linspace(0,3,N);

[T,U] = meshgrid(t,u);
DT = ones(size(T));
DU = f(T, U);

quiver(T,U,DT,DU)
xlim([0,5]); ylim([0,3])
hold on

[t, u] = ode45(f, [a, b], u0);

%
%% 4.

u0 = [0.1; 0.1];

f1 = @(u1, u2) 2*u1 - u1.*u2;               % We have two differential functions of (u1, u2)
f2 = @(u1, u2) u2 + 0.4*u1.*u2 - u2.^2;

f = @(t, u) [f1(u(1), u(2));
             f2(u(1), u(2))] % We create a system of equations to feed to ode45. It must be a function f(t,u), even if t is not involved (such as in this case).
[t, u] = ode45(f, [0, 8], u0)
close all
plot(t, u(:,1), t, u(:,2))
grid on

figure
U1 = u(:,1);
U2 = u(:,2);
plot(U1, U2)
hold on

u1 = linspace(-1, 4, 25);
u2 = linspace(-1, 3, 25);

[U1, U2] = meshgrid(u1, u2);

DU1 = f1(U1, U2);
DU2 = f2(U1, U2);

quiver(U1,U2,DU1, DU2, 1.5)
axis([-1, 4, -1, 3])

%
%% 4.

close all

f1 = @(u1, u2) cos(1 + u1 - u2.^2) + 0.3;
f2 = @(u1, u2) sin(0.3 + u1 + u2) - 0.4.*u1;

u0 = [1.5; 3];

f = @(t, u) [f1(u(1), u(2));
             f2(u(1), u(2))];

[t, u] = ode45(f, [0, 10], u0);
plot(t, u(:,1), t, u(:,2))
figure
plot(u(:,1), u(:,2))
hold on

u1 = linspace(0, 3, 25);
u2 = linspace(0, 3, 25);

[U1, U2] = meshgrid(u1, u2);
DU1 = f1(U1, U2);
DU2 = f2(U1, U2);

quiver(U1, U2, DU1, DU2)
axis([0,3,0,3])

%
%% 5.

% u   = u1
% u'  = u2
% u'' = u3

g = 9.81;
L = 0.1;
m = 0.1;
c = 0.2;

uprim = @(u1, u2) u2;
uprimprim = @(u1, u2) -g/L .*sin(u1) -c/m.*u2;

f = @(t, u) [uprim(u(1), u(2)); uprimprim(u(1), u(2))];
u0 = [pi/2; 0];
close all
[t, u] = ode45(f, [0, 3], u0);

u1 = u(:,1);
u2 = u(:,2);
plot(u1, u2, 'k')
hold on
scatter(u1, u2, 20, linspace(0, 3, length(t)), "filled")

xlabel('\phi')
ylabel('\omega')
grid on
axis([-2, 2, -15, 15])